/* 3D Navier-Stokes Solver <-> (Level Set <-> Fast March) -> Marching Cubes */
#include <iostream>
#include "typen.h"
#include "matrix.h"
#include "navier.h"
#include "convective.h"
#include <math.h>
#include <errno.h>
#include <stdlib.h>
#include <cstring>

using namespace std;

Navier::Navier(int argc,char* argv[]) 
{  
   DefinirCenaInicial(argc,argv);//cria as estruturas de dados e define os flags de obst�culo, inflow, fluido,..   
   /*kabs[i]=new double[gridp[i]+3]; typen.cpp 265 (Scene::Init())
   absolute grid positions, kabs =  -0.125  0  0.125  0.25  0.375  0.5  0.625  0.75  0.875  1  1.125  
                                    -0.125  0  0.125  0.25  0.375  0.5  0.625  0.75  0.875  1  1.125  
                                    -0.125  0  0.125  0.25  0.375  0.5  0.625  0.75  0.875  1  1.125  */
   S.DumpInfo(stdout);//Apresenta algumas informa��es a respeito da cena e dos par�metros inicais
}

void Navier::DefinirCenaInicial(int argc,char* argv[]) 
{ 
   //cout<<"\n Entrei na Definir Cena Inicial \n";
   if(S.alphatg==-1) 
      S.alphatg=S.alpha;//alpha:(doadora)donor-cell upwind; parameter for velocity in momentum equations
   S.Init();// initialize some precalculated values: e.g. kabs absolute grid positions, d distancia entre 
   //os pontos, cell widths; in typen.cpp line 228     
   //inicia os vetores que cont�m o tamanho das matrizes velocidade (U,V,W) e algumas outras tempor�rias
   // gridp[]: number of gridpoints in each co-ordinate direction cavity: 32 x 32 x 32
   //{0,33,0,33,0,33}; dim2[6]={0,33, 0,33, 0,33}; S.gridp = {NX, NY, NZ}   
   int dim[6]  = {1,S.gridp[0],    1,  S.gridp[1],    1,  S.gridp[2]}; //para criar matrizes NX x NY x NZ nehuma c�lula fantasma; usado para eq de Poisoon e para o Gradiente Biconjugado
   int dim2[6] = {0,S.gridp[0]+1,  0,  S.gridp[1]+1,  0,  S.gridp[2]+1}; //para criar matrizes (NX+2) x (NY+2) x (NZ+2) uma c�lula fantasma de cada lado
   int dim3[6] = {0, S.gridp[0]+2, 0,  S.gridp[1]+2,  0,  S.gridp[2]+2}; //para criar matrizes (NX+3) x (Ny+3) x (NZ+3) tr�s c�lulas fantasmas for convective terms 1 or 2, 
                                                      //deveria ser S.GC e n�o 2; usado par U,V,W, T ...
   //na fronteira peri�dica aumenta-se uma unidade para cada lado, se S.GC for 2
   U.Init(dim3);//cria a matriz 3D e inicia seus dados com 0
   V.Init(dim3);
   W.Init(dim3);
   flag.Init(dim3); //flagAnterior.Init(dim3);
   P.Init(dim2);
   F.Init(dim2);//some temporarly used arrays for velocity, temperature and level-set
   G.Init(dim2);
   H.Init(dim2);   
   RHS.Init(dim);//RHS is the Right Hand Side array of the poisson equation, pagina 20
   //do guia do usu�rio: RHS � o lado direito de laplaciano de p^n+1 = (1/delta_t) * grad . u^til  
   Phi.Init(dim3);   
   // arrays for BiCGSTAB and BiCGSTAB+JPC Solvers para a equa��o de poisson
   //if (S.Solver == 5 || S.Solver == 6)//parece que s� foi implementado este
   { 	
      int i,j,k;
      BiCG_r0.Init(dim); BiCG_rj.Init(dim); 
      BiCG_pj.Init(dim); BiCG_vj.Init(dim);
      BiCG_sj.Init(dim); BiCG_tj.Init(dim); 
      MatBuf.Init(dim2); KoeffPreCond.Init(dim);
      IJKALLLOOP flag[i][j][k] = AIR;
      Init_Velocity_Inflow();
      Create_slip_list(slip_list_size);
      Create_outflow_list(outflow_list_size);
      //slip_list_size =(int)slip_list.size();
      //cout<<"\n sliplist_size ="<<slip_list_size<<"\n";
      //outflow_list_size = (int)outflow_list.size();
      Prepare_Structure_for_Slip_Outflow_Cells(slip_list,slip_list_size,1,SLIP_ZERO_TWO_FREE_FACES);
      //cout<<"\n sliplist_size ="<<slip_list_size<<"\n";
      Prepare_Structure_for_Slip_Outflow_Cells(outflow_list,outflow_list_size,1,OUTFLOW_ZERO_TWO_FREE_FACES);
      Call_Complete_Cell();
      //cout<<"\n sliplist_size ="<<slip_list_size<<"\n";
      //for(int i_temp=0;i_temp<slip_list_size;i_temp++)
      //cout<<"\n Slip_list"<<slip_list[i_temp]<<"\n";
   }
}
void Navier::Create_slip_list(int &count)
{
   //cout<<"\n Entrei na slip_list \n";
   count = 0;
   for(int i=0;i<=S.gridp[0]+1;i++)
   {
      for(int j=0;j<=S.gridp[1]+1;j++)
      {
         for(int k=0;k<=S.gridp[2]+1;k++)
         {
            if((i==0) || (j==0) || (k==0) || (i==S.gridp[0]) || (k==S.gridp[2]))
            //if((i-j) >= 10)
            {
               if(flag[i][j][k] == AIR)
               {
                  flag[i][j][k] = SLIP;
                  int indice = i * (S.gridp[1]+2) * (S.gridp[2]+2) + j * (S.gridp[2]+2) + k;
                  slip_list[count]=indice;
                  count++;
               }
            }
         }
      }
   }
}
void Navier::Create_outflow_list(int &count)
{
    //cout<<"\n Entrei na outflow_list \n";
   count=0;
   for(int i=0;i<=S.gridp[0]+1;i++)
   {
      for(int j=0;j<=S.gridp[1]+1;j++)
      {
         for(int k=0;k<=S.gridp[2]+1;k++)
         {
            //if((j==0) && (((i - S.gridp[1])*(i - S.gridp[1])+
                                   //(k - 0.5 * S.gridp[2])*(k - 0.5 * S.gridp[2]))<16))
            if(i<0)
            {
               if(flag[i][j][k] == AIR)
               {   
                  flag[i][j][k] = OUTFLOW;
                  int indice = i * (S.gridp[1]+2) * (S.gridp[2]+2) + j * (S.gridp[2]+2) + k;
                  outflow_list[count] = indice;
                  count++;
               }
            }
         }
      }
   }
}
 
void Navier::Prepare_Structure_for_Slip_Outflow_Cells(/*std::vector <int>*/int *list, int &list_size, int list_offset,int list_two_free_faces)
{ //retiro das lista slip e outflow as c�lulas com mais de uma vizinha livre (air ou fluid), e defino o flag das 
  //c�luas retiradas como "list_two_free_faces", e os flgas das vizinhas livres como "flag+=list_offset"
   int d1[6] = {-1, 1, 0, 0, 0, 0};//deslocamentos indicando as c�lulas vizinhas nas 3 dire��es
   int d2[6] = {0, 0, -1, 1, 0, 0};
   int d3[6] = {0, 0, 0, 0, -1, 1};
   int i_after = 0;
   for(int i=0; i<list_size; i++)
   {
      int neighbors_free = 0;
      int indice = list[i];     
      int k_sl_out = indice % (S.gridp[2]+2);
      int j_sl_out = indice /(S.gridp[2]+2) % (S.gridp[1]+2);
      int i_sl_out = indice / ((S.gridp[1]+2)*(S.gridp[2]+2)); 
      int i_sl_out_1,j_sl_out_1,k_sl_out_1;
      int n1;
      for(int n=0; n<6; n++)
      {
         if((i_sl_out + d1[n]> 0) && (i_sl_out + d1[n]< S.gridp[0]+1) &&
            (j_sl_out + d2[n]> 0) && (j_sl_out + d2[n]< S.gridp[1]+1) &&
            (k_sl_out + d3[n]> 0) && (k_sl_out + d3[n]< S.gridp[2]+1))
         {
            if(flag[i_sl_out + d1[n]][j_sl_out + d2[n]][k_sl_out + d3[n]] == AIR)
            {
               neighbors_free++;
               if(neighbors_free == 1)
               {
                  i_sl_out_1 = i_sl_out + d1[n];
                  j_sl_out_1 = j_sl_out + d2[n];
                  k_sl_out_1 = k_sl_out + d3[n];
                  n1 = n;
               }
               else
               {
                  int temp_offset = flag[i_sl_out + d1[n]][j_sl_out + d2[n]][k_sl_out + d3[n]]%
                                    OFFSET_LIMIT;
                  if(temp_offset==0)
                     flag[i_sl_out + d1[n]][j_sl_out + d2[n]][k_sl_out + d3[n]] += list_offset;
                  if(neighbors_free == 2)
                  {
                     temp_offset = flag[i_sl_out_1][j_sl_out_1][k_sl_out_1]%
                                    OFFSET_LIMIT;
                     if(temp_offset==0)
                        flag[i_sl_out_1][j_sl_out_1][k_sl_out_1] += list_offset;
                     flag[i_sl_out][j_sl_out][k_sl_out] = list_two_free_faces;
                  }
               }
            }
         }
      }
      if(neighbors_free == 1)
      {
         list[i_after] = list[i];
         flag[i_sl_out][j_sl_out][k_sl_out] +=n1;
         //cout<<"\n list[i_after],slip="<<list[i_after]<<","<<slip_list[i_after]<<"\n";
         //cout<<"\n list[i],iafter="<<list[i]<<","<<i_after<<"\n";
         //cout<<"\n i,j,k ="<<i_sl_out<<","<<j_sl_out<<","<<k_sl_out<<"\n";
         //cout<<"\n n1="<<n1<<"\n";
         i_after ++;
      }
   }
   list_size = i_after;
   //cout<<"\n iafter="<<i_after<<"\n";
}

void Navier::Init_Velocity_Inflow()
{  //estou calculando Phi (inicial) aqui novamente (j� calculei em Container::SurfaceInicial, para n�o ter comunicar essa classe
   // com a classe LevelSet. Observe que quando alteramos Phi aqui, precisamos alterar lset.Phi em SurfaceInicial (Container)
   for(int i=0;i<=S.gridp[0]+1;i++)
   {
      float ic = i +0.5;
      for(int j=0;j<=S.gridp[1]+1;j++)
      {
         float jc = j +0.5;
         for(int k=0;k<=S.gridp[2]+1;k++)
         {
            float kc = k +0.5;
            double  signdisthor= sqrt((jc-16)*(jc-16) + (kc-16)*(kc-16))-5.0; 
            double  disthor = fmax(signdisthor,0);  
          if(i>0)Phi[i][j][k] = 1.0/S.gridp[0] * (sqrt(disthor*disthor+(ic-1.0)*(ic-1.0)))-1.0E-5 ; 
          else Phi[i][j][k] = 1.0/S.gridp[0] * fmax(-0.5,signdisthor);
	  //if(i<=S.gridp[0])Phi[i][j][k] = 1.0/S.gridp[0] * (sqrt(disthor*disthor+(i-S.gridp[0])*(i-S.gridp[0])))-1.0E-5 ; 
          //else Phi[i][j][k] = 1.0/S.gridp[0] * fmax(-1.0,signdisthor);
         }
      }
   }
   //int count=0;
   //Definindo as condicoes de fronteira iniciais:
   for (int j=0;j<=S.gridp[1]+1;j++)
   {
      for (int k=0;k<=S.gridp[2]+1;k++)
      {
         if(Phi[0][j][k] <= 0)
            {
               //U[S.gridp[0]+1][j][k] = 0.3; 
               //U[S.gridp[0]][j][k] = 0.3;
               //flag[S.gridp[0]+1][j][k] = INFLOW;
               U[0][j][k] = 0.3; 
               U[1][j][k] = 0.3;
               flag[0][j][k] = INFLOW;
               //int indice =  (S.gridp[2]+2) + j * (S.gridp[2]+2) + k;
               //   inflow_list[count] = indice;
               //   count++;
            }
      }
   } 
   //inflow_list_size == count;
}
void  Navier::Call_Complete_Cell() 
{
   int countBorderList = 0;
   int sizeVectorBorder = S.gridp[0] * S.gridp[1] * S.gridp[2] / 2;//verificar se este � o tamanho ideal
   int i_border_fluid[3][sizeVectorBorder],j_border_fluid[3][sizeVectorBorder];
   int k_border_fluid[3][sizeVectorBorder], par_border_fluid[3][sizeVectorBorder];
                                                //[0][] armazena as c�lulas FLUID com uma �nica face livre
                                                //[1][] armazena as c�lulas FLUID com 2 faces livres
                                                //[2][] armazena as c�lulas FLUID com 3 faces livres
   int index_border_list[3]; 
   //int sgridp12 = S.gridp[1] * S.gridp[2];//apenas para evitar calcular esse produto v�rias vezes
   for(int i=0; i<3; i++)
      index_border_list[i]=0; 
   for(int i=1;i<=S.gridp[0];i++)
   {
      for(int j=1;j<=S.gridp[1];j++)
      {
         for(int k=1;k<=S.gridp[2];k++) 
         {
             IFFLUID(flag[i][j][k])//criando a lista de c�lulas FLUID que s�o vizinhas (anterior) a uma c�lula AIR
            {  //observe que ainda dentro deste loop, classifico as c�lulas como AIR_FLUID_END ou AIR_TO_FLUID, logo a 
               //c�lula i-1 j� pode ter flag == AIR_FLUID_END
              
               int m=-1; 
               int n=0;
               if((i<=S.gridp[0]) && (ISAIR_OUTFLOW(flag[i+1][j][k])))//Porque preciso completar as c�lulas AIR_FLUID_END???
               //se entra aqui � porque a c�lula (i-1) � FLUID e a c�lula i � AIR, observe que a c�lula a ser inserida �
               //a c�lula (i-1) (FLUID)
                  {m ++; n=4;}//representa o n�mero (n-1) de faces livres (que presisam ter suas respectivas velocidades determinadas)
               if((j<=S.gridp[1]) && (ISAIR_OUTFLOW(flag[i][j+1][k])))
                  {m ++; n +=2;}
               if((j<=S.gridp[2]) && (ISAIR_OUTFLOW(flag[i][j][k+1]))) 
                  {m++; n++;}
               n=7-n;
               if(m>=0)
               {
                  i_border_fluid[m][index_border_list[m]]=i;
                  j_border_fluid[m][index_border_list[m]]=j;
                  k_border_fluid[m][index_border_list[m]]=k;
                  par_border_fluid[m][index_border_list[m]]=n;
                  index_border_list[m]++;
                  countBorderList++;
               }
            }
         }
      }
   }
   for(int p=0; p<3; p++)   
      for(int q=0; q<index_border_list[p]; q++)         
         Complete_Cell(i_border_fluid[p][q],j_border_fluid[p][q],k_border_fluid[p][q],par_border_fluid[p][q]);
}
void Navier::DoIt(Velocity gridCellUVW)
{  
   
   //cout<<"\n Entrei na DoIt \n";
   
   num_iter_Navier ++;
   bool flagsixfluidneighbors = false;
   bool Phiant,Phinew;
   bool Roant,Ronew;
   double res=0;
   int i,j,k;     
   S.delt=0.0;
   S.timestepmethod=0; 
   int countAirToFluid = 0;
   int countFluidToAir = 0;
   int countAirBorder = 0;
   for(int i=1;i<=S.gridp[0];i++)
   {
      for(int j=1;j<=S.gridp[1];j++)
      {
         /*if(Phi[i][j][1] > 0) 
         {
            Phiant = true;
            Roant = (Phi[i][j+1][1] <= 0) || (Phi[i+1][j][1] <= 0) || (Phi[i+1][j+1][1] <= 0);
         }
         else
         {
            Roant = true;
            Phiant = (Phi[i][j+1][1] > 0) || (Phi[i+1][j][1] > 0) || (Phi[i+1][j+1][1] > 0);
         }*/
         for(int k=1;k<=S.gridp[2];k++)     
         {   
            if(Phi[i][j][k] > 0)
            {
            //AIR          = 0x1000,     // =4096;   
            //AIR_TO_FLUID = 0x1F00,     // =7936;
            //AIR_FLUID_END = 0x1FF0,     // =8176; 
            //FLUID        = 0x1FFF,     // =8192;
            //ainda dentro dos for
            /*if(Phi[i][j][k+1] > 0) 
            {
               Phinew = true;
               Ronew = (Phi[i][j+1][k+1] <= 0) || (Phi[i+1][j][k+1] <= 0) || (Phi[i+1][j+1][k+1] <= 0);
            }
            else
            {
               Ronew = true;
               Phinew = (Phi[i][j+1][k+1] > 0) || (Phi[i+1][j][k+1] > 0) || (Phi[i+1][j+1][k+1] > 0);
            }
            if( Phiant ||  Phinew)
            {//o if acima funciona como uma nova verificacao de flag, ou seja, aqui eh o caso de n�o ser fluido*/
               IFFLUID(flag[i][j][k])
               {
                  countFluidToAir++;
                  int data = i*S.gridp[1]*S.gridp[2] + j*S.gridp[2] + k;  
                  //cout<<"\n Vou acessar a FluidToAirlist \n";
                  listFluidToAir.push_back(data);
               }
               /* if(Roant ||  Ronew) 
               {
                  countAirBorder++;
                  int data = i*S.gridp[1]*S.gridp[2] + j*S.gridp[2] + k;
                  //cout<<"\n Vou acessar a Airtoborderlist \n";
                  listAirBorder.push_back(data);
               }*/

            }
            else //eh fluido na nova verificacao
            {
               IFAIR(flag[i][j][k])//Se era AIR e passou a ser FLUID; estamos considerando que "sempre" existe uma 
               {                 //vizinha FLUID, caso n� haja envio uma mensagem, ????mas deveria tratar???
                  countAirToFluid++;
                  int data = i*S.gridp[1]*S.gridp[2] + j*S.gridp[2] + k;  
                  //cout<<"\n Vou acessar a AirtoFluidlist \n";
                  listAirToFluid.push_back(data);
                  //cout<<"\n AIR TO FLUID ="<<i<<","<<j<<","<<k<<"\n";
               }      
            } 
            //Phiant = Phinew;
            //Roant = Ronew;
         }//aqui termina os for 1<= x,y,z <= N
      }
   }
   //*****************************************************************************************************
   double uvmax = 0.0;
   int imax,jmax,kmax,control;
   /*IJKLOOP
   {
      //double div = 0.0;
      if((flag[i][j][k] > AIR_FLUID_END) || (flag[i-1][j][k] > AIR_FLUID_END))
         if(fabs(U[i][j][k]) > fabs(uvmax)) {uvmax = U[i][j][k]; imax =i;jmax=j; kmax=k; control = 1;}
       if((flag[i][j][k] > AIR_FLUID_END) || (flag[i][j-1][k] > AIR_FLUID_END))      
         if(fabs(V[i][j][k]) > fabs(uvmax)) {uvmax = V[i][j][k];  imax =i;jmax=j; kmax=k; control =2;}
       if((flag[i][j][k] > AIR_FLUID_END) || (flag[i][j][k-1] > AIR_FLUID_END))
         if(fabs(W[i][j][k]) > fabs(uvmax)) {uvmax = W[i][j][k];  imax =i;jmax=j; kmax=k; control =3;}
   }
   cout<<" uvmaxantes,i,j,k,ct ="<<uvmax<<","<<imax<<","<<jmax<<","<<kmax<<","<<control<<"\n";*/
   //*****************************************************************************************************
   int sizeList = (int)listAirToFluid.size();//n�mero de novas c�lulas fluido
   int *DirectList = new int[sizeList];
   int size_updated_List = 0;
   Prepare_Air_To_Fluid(DirectList,&size_updated_List); 
   //******************************************************************************************************
   /*double uvmaxproj = 0.0;
   for(int m=0; m<sizeList;m++)
   {
      int indice = reverseList[m];
      k = indice % S.gridp[2];
      j = (indice / S.gridp[2]) % S.gridp[1];
      i = indice / (S.gridp[1]*S.gridp[2]); 
      
         if(fabs(U[i][j][k]) > fabs(uvmaxproj)) {uvmaxproj = U[i][j][k]; imax =i;jmax=j; kmax=k; control = 1;}
         if(fabs(V[i][j][k]) > fabs(uvmaxproj)) {uvmaxproj = V[i][j][k];  imax =i;jmax=j; kmax=k; control =2;}
         if(fabs(W[i][j][k]) > fabs(uvmaxproj)) {uvmaxproj = W[i][j][k];  imax =i;jmax=j; kmax=k; control =3;}
   }
   if( uvmaxproj > uvmax)
      cout<<" ############################Inicializa��o aumentou valor############################### \n";
   cout<<" uvmaxantes,uvmaxdepois,i,j,k,ct ="<<uvmax<<","<<uvmaxproj<<","<<imax<<","<<jmax<<","<<
               kmax<<","<<control<<"\n";*/
   //*************************************************************************************************
      //percorrer a lista FinalList chamando a rotina NewFluidCellInitializer para cada elemento   
   //double divMaxAfterNewFluidInitializer = 0.0;
   //double div_1; 
   //cout<<"\n countFluidToAir="<<countFluidToAir<<"\n";   
   //cout << "\n Divergente maximo da lista AirToFluid apos NewFluidInitializer = " << divMaxAfterNewFluidInitializer;
   //substituir o loop abaixo por uma varrdeura das listas FluidToAir e AirToFluid e atualizar os flag's   
   vector<int>::iterator it = listFluidToAir.begin();
   for(int m=0; m<countFluidToAir; m++) 
   {
      int indice = *it;     
      k = indice % S.gridp[2];
      j = (indice / S.gridp[2]) % S.gridp[1];
      i = indice / (S.gridp[1]*S.gridp[2]);       
      flag[i][j][k] = AIR+flag[i][j][k]%OFFSET_LIMIT;
      P[i][j][k] = 0.0;
      it++;
   }
   //cout << "\n Vou deletar listFluidToAir";
     listFluidToAir.erase(listFluidToAir.begin(),listFluidToAir.end());
   //cout << "\n deletei listFluidToAir";
   /*if(flagsixfluidneighbors)
   {   
      //cout<<"\n Entrei flagsix \n";
      CompRHS();// Adams-Bashfort or Runge-Kutta Compute right hand side for poisson equation	 
      res=CompPoisson(starttime); // Solve poisson equation    
      AdapUVW(); // Compute new velocity (velocity correction by pressure-gradient)         
   }*/
 
   //##############################   MAIN BLOCK   ######################################
   S.delt=TimeStep(); // Compute new time step de acordo com a condi��o de CFL (ligeiramente modificada, um pouco mais restritiva)   
   //cout<<"\n sliplist_size Call_Update ="<<slip_list_size<<"\n";
   //for(int i_temp=0;i_temp<slip_list_size;i_temp++)
   //   cout<<"\n Slip_list"<<slip_list[i_temp]<<"\n";
   Call_Update_Slip_Outflow_Cells();
   
   //cout<<"\n Chamei CompFGH \n";
   CompFGH(gridCellUVW);//calculate temporary velocity fields (in space) due to chorin projection method (see User's Guide)   
   //cout<<"\n Chamei CompTUVW \n";
   CompTUVWfromFGH(S.timestepmethod);// Compute predicted velocities with 2nd order 
   //cout<<"\n Chamei RHS \n";
   //cout<<"\n countAirToFluid"<<countAirToFluid<<"\n";
   for(int m=0; m<size_updated_List; m++)
   {
      //int indice = finalList[m];
      int indice = DirectList[m];
      k = indice % S.gridp[2];
      j = (indice / S.gridp[2]) % S.gridp[1];
      i = indice / (S.gridp[1]*S.gridp[2]);       
      flag[i][j][k] = FLUID + flag[i][j][k]%OFFSET_LIMIT;
      //cout<<"\n i,j,k,flag="<<i<<","<<j<<","<<k<<","<< flag[i][j][k]<<"\n";
   }
   
   //listAirBorder.erase(listAirBorder.begin(),listAirBorder.end());
   //cout << "\n deletei listAirtoFluid";
   CompRHS();// Adams-Bashfort or Runge-Kutta Compute right hand side for poisson equation	 
   IJKLOOP
   {  if(ISAIR(flag[i][j][k]))
         P[i][j][k] = 0.0;
   } 
   //cout<<"\n Chamei Poisson \n"; 
   res=CompPoisson(/*starttime*/); // Solve poisson equation  	   
   /*double pmax = -1000; 
   double pmin = 1000;
   int imin, jmin, kmin /*,imax, jmax, kmax;
   for (int i=1;i<=S.gridp[0];i++)
      for (int j=1;j<=S.gridp[1];j++) 
         for (int k=1;k<=S.gridp[2];k++) 
         {  
            if(P[i][j][k]<pmin)
            {
               pmin = P[i][j][k];
               imin = i; jmin = j; kmin = k;
            }
           if(P[i][j][k]>pmax)
           {
              pmax = P[i][j][k];
              imax = i; jmax = j; kmax = k;
           }
         }*/
   //std::cout << "\t pmin = " << pmin << "("<< imin << ", "<< jmin << "," << kmin << ")";
   //std::cout << "\t pmax = " << pmax << "("<< imax << ", "<< jmax << "," << kmax << ")";
   //cout<<"\n Chamei AdapUVW \n"; 
   AdapUVW(); // Compute new velocity (velocity correction by pressure-gradient) 
  // cout<<"\n Sai AdapUVW \n";      
	S.delt_old=S.delt; // old time step needed for adams-bashfort          
	S.timestepmethod|=1;//used to mark the first time step for AB discretization and data output 
	S.timesteps++; //tag which defines if F,G,H arrays are from actual or previous time step     
}
//Considerando que j� foi criada uma lista com as c�lulas AirToFluid, criar uma lista para inicializa��o das novas 
//c�lulas fluido mantendo a divergencia nula
void Navier::Prepare_Air_To_Fluid(int *DirectList, int *size_updated_List)
{   
   vector<int>::iterator it;
   int i, j, k, indice; 
   unsigned int flagLevel;    
   int posicao, reverseIndex, countLevel, initFlagLevel[4], reverseIndexBound ;   
   int d1[6] = {-1, 1, 0, 0, 0, 0};//deslocamentos indicando as c�lulas vizinhas nas 3 dire��es
   int d2[6] = {0, 0, -1, 1, 0, 0};
   int d3[6] = {0, 0, 0, 0, -1, 1};
   countLevel = 0;//n�mero de n�veis a partir das c�lulas que continuam Air
   //n�vel 1 = novas c�lulas fluido que tem uma face em comum com uma que continua Air
   //n�vel i = novas c�lulas fluido que tem uma face em comum com outra dessas c�lulas do n�vel i-1 e n�o s�o adjacentes
   //          a c�lulas de n�veis mais baixos
   //se o deslocamento � limitado a 1 c�lula de uma itera��o para outra ent�o o n�vel <= 3 ???
   //Iremos come�ar a computar a velocidade nas c�lulas AirToFluid (divergencia=0) a partir das c�lulas de maior n�vel
   //(mais distantes de Air)
   flagLevel = AIR_FLUID_BEGIN;//n�vel correntemente tratado   
   int sizeList = (int)listAirToFluid.size();//n�mero de novas c�lulas fluido
   int *reverseList = new int[sizeList];//lista de novas c�lulas fluido ordenadas de forma que a de n�vel mais alto fique na frente   
   //finalList = new int[sizeList];//� a lista final de novas c�lulas que � ordenada por n�vel e dentro de um mesmo
   //n�vel as c�lulas mais pr�ximas das que eram fluido na itera��o anterior aparecem antes. Assim, quando formos
   //tratar uma c�lula para tornar sua diverg�ncia nula, teremos sempre uma componente da velocidade ainda n�o
   //definida pelo tratamento de c�lulas resolvidas anteriormente (tem-se sempre 1 grau de liberdade)
   reverseIndex = sizeList - 1;//ind�ce da lista reversa onde as c�lulas mais distantes de Air aparecem primeiro (na frente).
   reverseIndexBound = sizeList - 1;
   initFlagLevel[0] = sizeList;//initFlagLevel[i] indica a posicao do primeiro elemento da lista reversa no n�vel i+1. M�ximo de 4 n�veis
   do
   {
      //cout << "\n countLevel = " << countLevel << "     intiFlagLevel = " << initFlagLevel[countLevel];
      posicao = 1;//posi��o na lista original (listAirToFluid) 
      it = listAirToFluid.begin(); 
      for(int m=0; m<=reverseIndexBound; m++)//Porque n�o usar o pr�prio listAirToFluid.size()-1 ???
      {        
         indice = *it;         
         k = indice % S.gridp[2];
         j = (indice / S.gridp[2]) % S.gridp[1];
         i = indice / (S.gridp[1]*S.gridp[2]); 
         
         for(int n=0; n<6; n++)
         {  
            if(/*(i+d1[n] >=0) && (j+d2[n]>=0) && (k+d3[n]>=0) && 
               (i+d1[n] <=S.gridp[0]+1) && (j+d2[n]<=S.gridp[1]+1) && (k+d3[n]<=S.gridp[2]+1) && */
               (flag[i+d1[n]][j+d2[n]][k+d3[n]] <= flagLevel))//se houver uma vizinha em n�vel menor ou igual (na primeira itera�ao               
            {//do "do", � equivalente a verificar se existe alguma vizinha Air
               flag[i][j][k] = flagLevel + OFFSET_LIMIT; //passa ter o n�vel flagLevel+1 e ... 
               it--;//verificar a possibilidade de passar o it para a fun��o erase abaixo ???
               listAirToFluid.erase (listAirToFluid.begin()+posicao-1);// � retirado da lista AirToFluid e ...               
               reverseList[reverseIndex] = indice;//inserido na lista reverseList               
               reverseIndex--;//o �ndice da reverseList diminui porque a mesma � preenchida de tr�s para frente               
               posicao--;//pois quando um elemento � removido o pr�ximo elemento assume a posi��o dele, e uma vez
               
               break; //que estamos incrementando a posi��o a cada itera�a�o (logo abaixo) devemos compensar isto aqui
            }
         }
         it++;
         posicao++; 
         
      }
      
      countLevel++;      
      
      initFlagLevel[countLevel] = reverseIndex + 1;
      //cout << "\n countLevel = " << countLevel << "     intiFlagLevel = " << initFlagLevel[countLevel] << "\n\n\n\n\n";
      
      reverseIndexBound = reverseIndex; 
      flagLevel+=OFFSET_LIMIT; 
      if(flagLevel > AIR_FLUID_MIDDLE)
      {  
      cout<<"\n Unexpected number of levels among the new fluid cells\n";
      exit(0);
      }
      
   }while(listAirToFluid.size() > 0);
      
   //cout << "\n Nivel maximo das celula AirToFluid = " << countLevel;
   flagLevel = AIR_FLUID_END;
   for(int m=countLevel; m>0; m--)
   {        
      //cout << " m = " << m;
      int sizeListParcial = initFlagLevel[m-1] - initFlagLevel[m];
      orderAirToFluid(&(reverseList[initFlagLevel[m]]), sizeListParcial, &flagLevel,DirectList,size_updated_List);
      //cout << "\n Computei initFlagLevel[countLevel-1] - initFlagLevel[countLevel] = "  << initFlagLevel[countLevel-1] <<", "<< initFlagLevel[countLevel];
      /*int *a ;      
      a = orderAirToFluid(&(reverseList[initFlagLevel[m]]), sizeListParcial, &flagLevel);            
      for(int i = 0; i< sizeListParcial; i++)
      {
         finalList[initFlagLevel[m]+i] = a[i];             
      }//o vetor a est� sendo apresentado "corretamente" quando countLevel=1, o problema � quando countLevel � maior que 1
      */
   }        
}

//int* Navier::orderAirToFluid(int *list, int sizeList, unsigned int *flagLevel)
void Navier::orderAirToFluid(int *list, int sizeList, unsigned int *flagLevel,int *DirectList, int *size_updated_List)
{
   //cout<<"\n Entrei na orderAir \n";
   int i, j, k;
   int indice, countReverse = 0;
   int countDirect = *size_updated_List;   
   //int *directList = new int[sizeList];
   bool moved;
   //unsigned int flagLevel = AIR_FLUID_END;
   int d1[26] = {-1, 1, 0, 0, 0, 0,   -1,-1,-1,-1,  0,  0, 0, 0,  1, 1, 1, 1,   -1, -1, -1, -1,  1,  1, 1, 1};
   int d2[26] = {0, 0, -1, 1, 0, 0,   -1, 0, 0, 1, -1, -1, 1, 1, -1, 0, 0, 1,   -1, -1,  1,  1, -1, -1, 1, 1};
   int d3[26] = {0, 0, 0, 0, -1, 1,    0,-1, 1, 0, -1,  1,-1, 1,  0,-1, 1, 0,   -1,  1, -1,  1, -1,  1,-1, 1};   
                  // viz. p/faces         //viz. p/ arestas                                 //viz. p/ vertices
   int neighborhoodlevel = 0;
   int bound[4] = {0, 6, 18, 26};
   //double divmax = 0.0;
   do         
   {   
      /* cout<<"\n sizeList, neighborhoodlevel="<<sizeList<<","<<neighborhoodlevel<<"\n";
       for(int l=0; l<sizeList; l++)
      {
         indice = list[l];         
         k = indice % S.gridp[2];
         j = (indice / S.gridp[2]) % S.gridp[1];
         i = indice / (S.gridp[1]*S.gridp[2]); 
         cout<<"\n i,j,k ="<<i<<","<<j<<","<<k<<"\n";
         cout<<"\n flag(i,j,k)="<<flag[i][j][k]<<"\n";
      }*/
      for(int l=0; l<sizeList; l++)
      {
         indice = list[l];         
         k = indice % S.gridp[2];
         j = (indice / S.gridp[2]) % S.gridp[1];
         i = indice / (S.gridp[1]*S.gridp[2]);                 
         for(int n=bound[neighborhoodlevel]; n<bound[neighborhoodlevel+1]; n++)
         {
            //cout<<"\n id1,jd2,kd3,n ="<<i+d1[n]<<","<<j+d2[n]<<","<<k+d3[n]<<","<<n<<"\n";
            moved = false;
            if(/*(i+d1[n] >=0) && (j+d2[n]>=0) && (k+d3[n]>=0) && 
               (i+d1[n] <=S.gridp[0]+1) && (j+d2[n]<=S.gridp[1]+1) && (k+d3[n]<=S.gridp[2]+1) && */
               (flag[i+d1[n]][j+d2[n]][k+d3[n]] >= *flagLevel))
            {
               flag[i][j][k] = *flagLevel - OFFSET_LIMIT;
               //cout<<"\n i,j,k ="<<i<<","<<j<<","<<k<<"\n";
               NewFluidCellInitializer(i, j, k);
               //double div_1 = U[i+1][j][k] + V[i][j+1][k] + W[i][j][k+1] - (U[i][j][k] + V[i][j][k] + W[i][j][k]);
               //if(fabs(div_1) > fabs(divmax)) divmax = div_1;
               //cout<<"\n UVW[i,j,k] ="<<U[i][j][k]<<","<<V[i][j][k]<<","<<W[i][j][k]<<"\n";
               DirectList[countDirect] = indice;
               countDirect++;
               moved = true;
               break;               
            }            
         }         
         if(!moved)
         {
            list[countReverse] = list[l];
            countReverse++;
         }
      }   
      if (sizeList == countReverse)
      {
         if(neighborhoodlevel<2)
         {
            neighborhoodlevel++;
            countReverse = 0;
         }
         else
	 {
	    double epsilon = 1E-8;
            for(int itemp = 0;itemp < sizeList;itemp++)
            {
	       double  phimin = 0;
               indice = list[itemp];         
               k = indice % S.gridp[2];
               j = (indice / S.gridp[2]) % S.gridp[1];
               i = indice / (S.gridp[1]*S.gridp[2]); 
               cout<<"\n I,J,K ="<< i <<","<< j <<","<< k <<"\n";
               //cout<<"FLAG ="<<flag[i][j][k]<<"\n";
               for(int deltai=0; deltai<2;deltai++)
                  for(int deltaj=0; deltaj<2;deltaj++)
                     for(int deltak=0; deltak<2;deltak++)
                     {
                         //cout<<"\ni,j,k ="<<i+deltai<<","<<j+deltaj<<","<<k+deltak<<"\n";
                         //cout<<"flag ="<<flag[i+deltai][j+deltaj][k+deltak]<<"\n";
			 if(Phi[i+deltai][j+deltaj][k+deltak]<phimin) phimin= Phi[i+deltai][j+deltaj][k+deltak];
                     }
               //cout<<"\ni,j,k ="<<i<<","<<j<<","<<k<<"\n";
               //cout<<"neilevel,sizelist ="<<neighborhoodlevel<<","<<sizeList<<"\n";
	     for(int deltai=0; deltai<2;deltai++)
                  for(int deltaj=0; deltaj<2;deltaj++)
                     for(int deltak=0; deltak<2;deltak++)
                     {
                         //cout<<"\ni,j,k ="<<i+deltai<<","<<j+deltaj<<","<<k+deltak<<"\n";
                         //cout<<"flag ="<<flag[i+deltai][j+deltaj][k+deltak]<<"\n";
			 Phi[i+deltai][j+deltaj][k+deltak] -= (phimin -epsilon);
		     }
            }
            cout << "\n unexpected fluid bubbles";
            //exit(0);
         }
      }
      else
      {
         neighborhoodlevel = 0;
         sizeList=countReverse;
         countReverse = 0;
         *flagLevel-=OFFSET_LIMIT;
      }
      
   }while(sizeList>0);
   *size_updated_List=countDirect;
   //cout<<"\n Sai da orderAir \n";
  
   //cout<<"\n DIVMAXINITIALIZER ------------------------>"<<divmax<<"\n";
}


void Navier::CompF(int i, int j, int k)//calculate temporary U velocity
{
   int border = true;
   double a0,a11,a12, b0, b11, b12;//,c1,c2,d1,d2;
   double du2dx, duvdy, duwdz, d2udx2, d2udy2, d2udz2;
   double Gx;
   Gx = S.g[0]/S.froude;
   S.nu = 0.0;
   //std::cout << "\n ******S.nu*******" << S.nu << "\n";
   //std::cout << "\n Antes do erro 1 \n";
   if((i>=2) && (i<=S.gridp[0]-1) && ISFLUID(flag[(i-1)][j][k]) && ISNOTAIR(flag[(i+1)][j][k]) && ISNOTAIR(flag[(i-2)][j][k]))
   {
      a0=U[i-2][j][k];
      b0=U[i+2][j][k];
      border = false;
   }
   //std::cout << "\n Antes do erro 2 \n";
   du2dx= DUU(a0,U[i-1][j][k],U[i][j][k],U[i+1][j][k],b0,DX[i-1],DX[i],DX[i+1],DX[i+2],S.alpha,border);
   //
   border=false;
   if((j>=2) && (ISFLUID(flag[i][j-1][k]))&&(ISNOTAIR(flag[i][j-2][k])))
   {
      a0=U[i][j-2][k];
      a11=U[i][j-1][k];
   }
   else
   {
      border=true;
      if(ISNOTAIR(flag[i][j-1][k]))
      {
         a11=U[i][j-1][k];
      }
      else
      {
         a11=U[i][j][k];
      }
   }
   if((j<=S.gridp[1]-1) && (ISNOTAIR(flag[i][j+2][k])) && (ISFLUID(flag[i][j+1][k])))
   {
      b0=U[i][j+2][k];
      b11=U[i][j+1][k];
   }
   else
   {
      border= true;
      if(ISNOTAIR(flag[i][j+1][k]))
      {
         b11=U[i][j+1][k];
      }
      else
      {
         b11=U[i][j][k];
      }
   }
   duvdy=DUV(a0,a11,U[i][j][k],b11,b0,(DX[i-1]*V[i][j][k]+DX[i]*V[i-1][j][k])/(DX[i]+DX[i-1]),
         (DX[i-1]*V[i][j+1][k]+DX[i]*V[i-1][j+1][k])/(DX[i]+DX[i-1]),DY[j-2],DY[j-1],DY[j],DY[j+1],DY[j+2],S.alpha,border);
                 
   border=false;
   if((k>=2) && (ISNOTAIR(flag[i][j][k-2])) && (ISFLUID(flag[i][j][k-1])))
   {
      a0=U[i][j][k-2];
      a12=U[i][j][k-1];
   }
   else
   {
      border=true;
      if((ISNOTAIR(flag[i][j][k-1])))
      {
         a12=U[i][j][k-1];
      }
      else
      {
        a12=U[i][j][k];
      }
   }
   if((k<=S.gridp[2]-1) && (ISNOTAIR(flag[i][j][k+2])) && (ISFLUID(flag[i][j][k+1])))
   {
      b0=U[i][j][k+2];
      b12=U[i][j][k+1];
  }
  else
  {
     border=true;
     if((ISNOTAIR(flag[i][j][k+1])))
     {
        b12=U[i][j][k+1];
     }
     else
     {
        b12=U[i][j][k];
     }
  }
 
  duwdz=DUV(a0,a12,U[i][j][k],b12,b0,(DX[i-1]*W[i][j][k]+DX[i]*W[i-1][j][k])/(DX[i]+DX[i-1]),
  (DX[i-1]*W[i][j][k+1]+DX[i]*W[i-1][j][k+1])/(DX[i]+DX[i-1]),DZ[k-2],DZ[k-1],DZ[k],DZ[k+1],DZ[k+2],S.alpha,border);

 

  d2udx2=DD (a0,U[i][j][k],U[i+1][j][k],0,i); //??? a0
  d2udy2=DDS(a11,U[i][j][k],b11,1,j);
  d2udz2=DDS(a12,U[i][j][k],b12,2,k);
  F[i][j][k] = S.nu*(d2udx2+d2udy2+d2udz2) - du2dx-duvdy-duwdz + Gx;
  if(isnan(F[i][j][k]) || isinf(F[i][j][k]))
               {
                  cout<<"***CompF*****i,j,k="<<i<<","<<j<<","<<k<<"\n";
                  exit(0);
               }
  /*if(i==3 && j==10 && k==8)
  {
     cout<<"***********************3,10,8,du2dx,duvdy,duwdz ="<<  du2dx <<","<< duvdy<<","<< duwdz <<"\n";
     cout<<"***********************3,10,8,du2dx,duvdy,duwdz ="<<  d2udx2 <<","<< d2udy2<<","<< d2udz2 <<"\n";
  }*/
  /*for(int idel=-1;idel<=1;idel++) 
  {
    for(int jdel=-1;jdel<=1;jdel++)
    {
      for(int kdel=-1;kdel<=1;kdel++)
      {
        if(((ISNOTAIR(flag[i+idel][j+jdel][k+kdel])) ||
           (ISNOTAIR(flag[i+idel][j+jdel-1][k+kdel])))&&
           (fabs(U[i+idel][j+jdel][k+kdel]) > max_temp))
              max_temp = fabs(U[i+idel][j+jdel][k+kdel]);
      }
    }
  }
  if(temp > max_temp)
  {
    // cout<<"************************************************************************** \n";
    //cout<<"***CompF***i,j,k,temp,max_temp"<<i<<","<<j<<","<<k<<","<<temp<<","<<max_temp<<"\n";
  }*/
}

void Navier::CompG(int i, int j, int k)//calculate temporary V velocity
{
                  int border = true;
                  double a0,a11,a12, b0, b11, b12;
                  double dv2dy, duvdx, dvwdz, d2vdx2, d2vdy2, d2vdz2;
                  double  Gy;
                  Gy = S.g[1]/S.froude;
                  S.nu = 0.0;
                 //std::cout << "\n ******S.nu*******" << S.nu << "\n";
                 

                 if((j>=2)  && ISNOTAIR(flag[(i)][j-2][k]) && ISNOTAIR(flag[i][j+1][k]))
                 {
                    border = false;
                    a0=V[i][j-2][k];
                    b0=V[i][j+2][k];
      
                 }
                 dv2dy= DUU(a0,V[i][j-1][k],V[i][j][k],V[i][j+1][k],b0,DY[j-1],DY[j],DY[j+1],DY[j+2],S.alpha,border);
                 
                border=false;
                if((i>=2) && (ISNOTAIR(flag[i-2][j][k])) && (ISFLUID(flag[i-1][j][k])))
                {
                  a0=V[i-2][j][k];
                  a11=V[i-1][j][k];
                }
                else
                {
                  border=true; 
                  if((ISNOTAIR(flag[i-1][j][k])))
                  {
                     a11=V[i-1][j][k];
                  }
                  else
                  {
                     a11=V[i][j][k];
                  }
                }
                if((i<=S.gridp[0]-1) && (ISNOTAIR(flag[i+2][j][k])) && (ISFLUID(flag[i+1][j][k])))
                {
                  b0=V[i+2][j][k];
                  b11=V[i+1][j][k];
                }
                else
                {
                  border=true;
                  if((ISNOTAIR(flag[i+1][j][k])))
                  {
                    b11=V[i+1][j][k];
                  }
                  else
                  {
                    b11=V[i][j][k];
                  }
                }
                
                
                duvdx=DUV(a0,a11,V[i][j][k],b11,b0,(DY[j-1]*U[i][j][k]+DY[j]*U[i][j][k])/(DY[j]+DY[j-1]),
                (DY[j-1]*U[i+1][j][k]+DY[j]*U[i+1][j][k])/(DY[j]+DY[j-1]),DX[i-2],DX[i-1],DX[i],DX[i+1],DX[i+2],S.alpha,border);
                

                border=false;
                if((k>=2) && (ISNOTAIR(flag[i][j][k-2])) && (ISFLUID(flag[i][j][k-1])))
                {
                  a0=V[i][j][k-2];
                  a12=V[i][j][k-1];
                }
                else
                {
                   border=true;
                   if((ISNOTAIR(flag[i][j][k-1])))
                   {
                      a12=V[i][j][k-1];
                   }
                   else
                   {
                      a12=V[i][j][k];
                   }
                }
                if((k<=S.gridp[2]-1) && (ISNOTAIR(flag[i][j][k+2])) && (ISFLUID(flag[i][j][k+1])))
                {
                   b0=V[i][j][k+2];
                   b12=V[i][j][k+1];
                }
                else
                {
                   border=true;
                   if((ISNOTAIR(flag[i][j][k+1])))
                   {
                      b12=V[i][j][k+1];
                   }
                   else
                   {
                      b12=V[i][j][k];
                   }
                }
                
                dvwdz=DUV(a0,a12,V[i][j][k],b12,b0,(DY[j-1]*W[i][j][k]+DY[j]*W[i][j-1][k])/(DY[j]+DY[j-1]),
               (DY[j-1]*W[i][j][k+1]+DY[j]*W[i][j-1][k+1])/(DY[j]+DY[j-1]),DZ[k-2],DZ[k-1],DZ[k],DZ[k+1],DZ[k+2],S.alpha,border);

                 d2vdy2=DD (V[i][j-1][k],V[i][j][k],V[i][j+1][k],1,j);
                 d2vdx2=DDS(a11,V[i][j][k],b11,0,i);
                 d2vdz2=DDS(a12,V[i][j][k],b12,2,k);
		
                G[i][j][k] = S.nu*(d2vdx2+d2vdy2+d2vdz2) - dv2dy-duvdx-dvwdz + Gy;
                if(isnan(G[i][j][k]) || isinf(G[i][j][k]))
               {
                  cout<<"***CompG*****i+1,j,k="<<i<<","<<j<<","<<k<<"\n";
                  exit(0);
               }
                 double max_temp = 0;
                 for(int idel=-1;idel<=1;idel++) 
                 {
                     for(int jdel=-1;jdel<=1;jdel++)
                     {
                        for(int kdel=-1;kdel<=1;kdel++)
                        {
                          if(((ISNOTAIR(flag[i+idel][j+jdel][k+kdel])) ||
                             (ISNOTAIR(flag[i+idel][j+jdel-1][k+kdel])))&&
                             (fabs(V[i+idel][j+jdel][k+kdel]) > max_temp))
                                max_temp = fabs(V[i+idel][j+jdel][k+kdel]);
                       }
                    }
                  }
                  
                     //cout<<"************************************************************************** \n";
                     //cout<<"***CompG***i,j+1,k,max_temp"<<i<<","<<j<<","<<k<<","<<temp<<","<<max_temp<<"\n";
                 
               
}

void Navier::CompH(int i, int j, int k)//calculate temporary W velocity
{
                  int border=false;
                  double a0,a11,a12, b0, b11, b12;
                  double dw2dz, duwdx, dvwdy, d2wdx2, d2wdy2, d2wdz2;
                  double Gz;
                  S.nu = 0.0;
                  Gz = S.g[2]/S.froude;
                 if((k>=2) && (k<=S.gridp[2]-1) && ISNOTAIR(flag[(i)][j][k-2]) && ISFLUID(flag[i][j][k+1]))
                 {
                    border=false;
                    a0=W[i][j][k-2];
                    b0=W[i][j][k+2];
                 }
                 dw2dz= DUU(a0,W[i][j][k-1],W[i][j][k],W[i][j][k+1],b0,DZ[k-1],DZ[k],DZ[k+1],DZ[k+2],S.alpha,border);
                
                border = false;
                if((i>=2) && (ISNOTAIR(flag[i-2][j][k])) && (ISFLUID(flag[i-1][j][k])))
                {
                   a0=W[i-2][j][k];
                   a11=W[i-1][j][k];
                }
                else
                {
                  border = true;
                  if((ISNOTAIR(flag[i-1][j][k])))
                  {
                     a11=W[i-1][j][k];
                  }
                  else
                  {
                     a11=W[i][j][k];
                  }
                }
                if((i<=S.gridp[0]-1) && (ISNOTAIR(flag[i+2][j][k])) && (ISFLUID(flag[i+1][j][k])))
                {
                  b0=W[i+2][j][k];
                  b11=W[i+1][j][k];
                }
                else
                {
                   border=true;
                   if((ISNOTAIR(flag[i+1][j][k])))
                   {
                      b11=W[i+1][j][k];
                   }
                   else
                   {
                      b11=W[i][j][k];
                   }
                }
                duwdx=DUV(a0,a11,W[i][j][k],b11,b0,(DZ[k-1]*U[i][j][k]+DZ[k]*U[i][j][k-1])/(DZ[k]+DZ[k-1]),
                (DZ[k-1]*U[i+1][j][k]+DZ[k]*U[i+1][j][k-1])/(DZ[k]+DZ[k-1]),DX[i-2],DX[i-1],DX[i],DX[i+1],DX[i+2],S.alpha,border);

                 border=false;
                 if((j>=2) && (ISNOTAIR(flag[i][j-2][k])) && (ISFLUID(flag[i][j-1][k])))
                 {
                    a0=W[i][j-2][k];
                    a12=W[i][j-1][k];
                 }
                 else
                 {
                    border=true;
                    if((ISNOTAIR(flag[i][j-1][k])))
                    {
                       a12=W[i][j-1][k];
                    }
                    else
                    {
                       a12=W[i][j][k];
                    }
                  }
                  if((j<=S.gridp[1]-1) && (ISNOTAIR(flag[i][j+2][k])) && (ISFLUID(flag[i][j+1][k])))
                  {
                    b0=W[i][j+2][k];
                    b12=W[i][j+1][k];
                  }
                  else
                  {
                    border=true;
                    if((ISNOTAIR(flag[i][j+1][k])))
                    {
                      b12=W[i][j+1][k];
                    }
                    else
                    {
                      b12=W[i][j][k];
                    }
                 }
                 
                 dvwdy=DUV(a0,a12,W[i][j][k],b12,b0,(DZ[k-1]*V[i][j][k]+DZ[k]*V[i][j][k-1])/(DZ[k]+DZ[k-1]), 
                 (DZ[k-1]*V[i][j+1][k]+DZ[k]*V[i][j+1][k-1])/(DZ[k]+DZ[k-1]),DY[j-2],DY[j-1],DY[j],DY[j+1],DY[j+2],S.alpha,border);


                 d2wdz2=DD (W[i][j][k-1],W[i][j][k],W[i][j][k+1],2,k);
                 d2wdx2=DDS(a11,W[i][j][k],b11,0,i);
                 d2wdy2=DDS(a12,W[i][j][k],b12,1,j);
		
                 H[i][j][k] = S.nu*(d2wdx2+d2wdy2+d2wdz2) - dw2dz - duwdx - dvwdy + Gz; 
                  if(isnan(H[i][j][k]) || isinf(H[i][j][k]))
                 {
                     cout<<"***CompH*****i,j,k="<<i<<","<<j<<","<<k<<"\n";
                     exit(0);
                 }
                 double max_temp = 0;
                 for(int idel=-1;idel<=1;idel++) 
                 {
                     for(int jdel=-1;jdel<=1;jdel++)
                     {
                        for(int kdel=-1;kdel<=1;kdel++)
                        {
                          if(((ISNOTAIR(flag[i+idel][j+jdel][k+kdel])) ||
                             (ISNOTAIR(flag[i+idel][j+jdel][k-1+kdel])))&&
                             (fabs(V[i+idel][j+jdel][k+kdel]) > max_temp))
                                max_temp = fabs(W[i+idel][j+jdel][k+kdel]);
                       }
                    }
                  }
                  
                     //cout<<"************************************************************************** \n";
                     //cout<<"***CompG***i,j+1,k,max_temp"<<i<<","<<j<<","<<k<<","<<temp<<","<<max_temp<<"\n";
                
}
void Navier::CompFM1(Velocity gridCellUVW, int i, int j, int k, int gama)//calculate temporary U velocity
{
              //cout<<"\n ****************** Entrei na CompFM1 ******************************** \n";
              //cout<<"\n *****************************"<<i<<","<< j<<","<<k<<"***********************\n";
	      int alpha,beta;
	      double e0,e1,e2,c1,c2,d1,d2,f1,f2,g1,g2;
	      double coef0,coef1,coef2;
         double d2udx2, d2udy2, d2udz2;
         double Gx;
         Gx = S.g[0]/S.froude;
         S.nu = 0.0;
         double h = 1.0/S.gridp[0];
         DX[0] =h; DY[0]=h; DZ[0]=h;
         DX[i] =h; DY[j]=h; DZ[k]=h;
         DX[i+1] =h; DY[j+1]=h; DZ[k+1]=h;
         DX[i-1] =h; DY[j-1]=h; DZ[k-1]=h;
	      alpha = 2*gama -1;
	      beta = alpha + gama;

	      coef0 = alpha*U[i+gama][j][k];
	      if(coef0>=0)
         {
            coef0 *=(-1.0/DX[i]);
            e0 = U[i+1-gama][j][k];
            coef1 = V[i][j+1][k]+V[i][j][k];
            coef2 = W[i][j][k+1]+W[i][j][k];
	      }
	     else
	     {
            coef0 *=(1.0/DX[i+alpha]);
            e0 = U[i+gama][j][k];
            /*if((i<S.gridp[0]) || (alpha == -1))
            {
               e0 =.25*(gridCellUVW.gridU(i+beta,j,k)+gridCellUVW.gridU(i+beta,j+1,k) + 
               gridCellUVW.gridU(i+beta,j,k+1)+gridCellUVW.gridU(i+beta,j+1,k+1));
            }
            else
            {
               e0 = U[i+1][j][k];//Se o fluxo está vindo de uma ghost cell ela só pode ser inflow 
            }*/

            if(ISNOTAIR_OUTFLOW(flag[i+alpha][j-1][k]))
            {
               c1 = V[i+alpha][j][k];
            }
            else 
               c1 = V[i][j][k];
            /*{
               if(i<S.gridp[0] || (alpha == -1))
               {	
                  c1 = .25*(gridCellUVW.gridV(i+gama,j,k)+gridCellUVW.gridV(i+beta,j,k) + 
                  gridCellUVW.gridV(i+gama,j,k+1)+gridCellUVW.gridV(i+beta,j,k+1));
               }
               else
               {
                  c1 = .5*(gridCellUVW.gridV(i+1,j,k)+ gridCellUVW.gridV(i+1,j,k+1));
               }
            }*/
			
            if(ISNOTAIR_OUTFLOW(flag[i+alpha][j+1][k]))
            {
               d1 = V[i+alpha][j+1][k];
            }
            else 
               d1 = V[i][j+1][k];
            /*{
               if(i<S.gridp[0] || (alpha == -1))
               {
                  d1 =.25*(gridCellUVW.gridV(i+gama,j+1,k)+gridCellUVW.gridV(i+beta,j+1,k) + 
                  gridCellUVW.gridV(i+gama,j+1,k+1)+gridCellUVW.gridV(i+beta,j+1,k+1));
               }
               else
               {
                  d1 = .5*(gridCellUVW.gridW(i+1,j+1,k)+ gridCellUVW.gridW(i+1,j+1,k+1));
               }
            }*/
            coef1 = c1+d1;
		
            if(ISNOTAIR_OUTFLOW(flag[i+alpha][j][k-1]))
            {
               c2 = W[i+alpha][j][k];
            }
            else 
               c2 = W[i][j][k];
            /*{
               if(i<S.gridp[0] || (alpha == -1))
               {	
                  c2 = .25*(gridCellUVW.gridW(i+gama,j,k)+gridCellUVW.gridW(i+beta,j,k) + 
                        gridCellUVW.gridW(i+gama,j+1,k)+gridCellUVW.gridW(i+beta,j+1,k));
               }
               else
               {
                  c2 = .5*(gridCellUVW.gridW(i+1,j,k)+ gridCellUVW.gridW(i+1,j+1,k));
               }
            }*/
			
            if(ISNOTAIR_OUTFLOW(flag[i+alpha][j][k+1]))
            {
               d2 = W[i+alpha][j][k+1];
            }
            else 
               d2 = W[i][j][k+1];
            /*{
               if(i<S.gridp[0] || (alpha == -1))
               {
                  d2 =.25*(gridCellUVW.gridW(i+gama,j,k+1)+gridCellUVW.gridW(i+beta,j,k+1) + 
                  gridCellUVW.gridW(i+gama,j+1,k+1)+gridCellUVW.gridW(i+beta,j+1,k+1));
               }
               else
               {
                  d2 = .5*(gridCellUVW.gridW(i+1,j,k+1)+ gridCellUVW.gridW(i+1,j+1,k+1));
               }
            }*/
		  coef2 = c2 +d2;

	      }
	      if((coef1)>= 0)
	      {
            coef1 *=(-1.0/(DY[j-1]+DY[j]));
            if(ISFLUID_INFLOW(flag[i][j-1][k]))
            {
               f1=e1=U[i+gama][j-1][k];
            }
            else
            {
               f1=e1=U[i+gama][j][k];
               /*f1=U[i+gama][j][k];
               e1 =.25*(gridCellUVW.gridU(i+gama,j-1,k)+gridCellUVW.gridU(i+gama,j-1,k+1) + 
                  gridCellUVW.gridU(i+gama,j,k)+gridCellUVW.gridU(i+gama,j,k+1));*/
            }
            if(ISFLUID_INFLOW(flag[i][j+1][k]))
            {
               g1=U[i+gama][j+1][k];
            }
            else
            {
               g1=U[i+gama][j][k];
            }
         }
         else
         {
            coef1 *=(1.0/(DY[j+1]+DY[j]));
            if(ISFLUID_INFLOW(flag[i][j+1][k]))
            {
               g1=e1=U[i+gama][j+1][k];
            }
            else
            {	
               g1=e1=U[i+gama][j][k];
               /*g1=U[i+gama][j][k];
               if(j<S.gridp[1])
               {
                  e1 =.25*(gridCellUVW.gridU(i+gama,j+1,k)+gridCellUVW.gridU(i+gama,j+1,k+1) + 
                           gridCellUVW.gridU(i+gama,j+2,k)+gridCellUVW.gridU(i+gama,j+2,k));
               }
               else
               {
               e1 =.5*(gridCellUVW.gridU(i+1,j+1,k)+gridCellUVW.gridU(i+1,j+1,k+1));
               }*/
            }
            if(ISFLUID_INFLOW(flag[i][j-1][k]))
            {  
               f1=U[i+gama][j-1][k];
            }
            else
            {
               f1=U[i+gama][j][k];
            }
	      }
	      if((coef2)>=0)
	      {
            coef2 *=(-1.0/(DZ[k-1]+DZ[k]));
            if(ISFLUID_INFLOW(flag[i][j][k-1]))
            {
               f2=e2=U[i+gama][j][k-1];
            }
            else
            {
               f2=e2=U[i+gama][j][k];
               /*f2=U[i+gama][j][k];
               e2 =.25*(gridCellUVW.gridU(i+gama,j,k-1)+gridCellUVW.gridU(i+gama,j+1,k-1) + 
                     gridCellUVW.gridU(i+gama,j,k)+gridCellUVW.gridU(i+gama,j+1,k));*/
            }
            if(ISFLUID_INFLOW(flag[i][j][k+1]))
            {
               g2=U[i+gama][j][k+1];
            }
            else
            {
		        g2=U[i+gama][j][k];
            }
         }
         else
         {
            coef2 *=(1.0/(DZ[k+1]+DZ[k]));
            if(ISFLUID_INFLOW(flag[i][j][k+1]))
            {
               g2=e2=U[i+gama][j][k+1];
            }
            else
            {	
               g2=e2=U[i+gama][j][k];
               /*g2=U[i+gama][j][k];
               if(k<S.gridp[2])
               {
                  e2 =.25*(gridCellUVW.gridU(i+gama,j,k+1)+gridCellUVW.gridU(i+gama,j+1,k+1) + 
                  gridCellUVW.gridU(i+gama,j,k+2)+gridCellUVW.gridU(i+gama,j+1,k+2));
               }
               else
               {
                  e2 =.5*(gridCellUVW.gridU(i+1,j,k+1)+gridCellUVW.gridU(i+1,j+1,k+1));
               }*/
            }
            if(ISFLUID_INFLOW(flag[i][j][k+1]))
            {
               f2=U[i+gama][j][k-1];
            }
            else
            {
               f2=U[i+gama][j][k];
            }
         }
         double derconv= coef0 *(e0-U[i+gama][j][k])+ coef1*(e1-U[i+gama][j][k]) + coef2*(e2-U[i+gama][j][k]);

         d2udx2=DD(U[i][j][k],U[i+gama][j][k],U[i+1][j][k],0,i);
         d2udy2=DDS(f1,U[i+gama][j][k],g1,1,j);
         d2udz2=DDS(f2,U[i+gama][j][k],g2,2,k);

         F[i+gama][j][k] = S.nu*(d2udx2+d2udy2+d2udz2) - derconv + Gx; 
         
         /*double coefsum = (fabs(coef0)+fabs(coef1)+fabs(coef2))*S.delt ;
         if(fabs(coef0) > (UVWmax/h)) cout<<"coef0 > UVWMax="<<coef0<<","<<UVWmax/h<<"\n";
         if(fabs(coef1) > (UVWmax/h)) cout<<"coef1 > UVWMax="<<coef1<<","<<UVWmax/h<<"\n";
         if(fabs(coef2) > (UVWmax/h)) cout<<"coef2 > UVWMax="<<coef2<<","<<UVWmax/h<<"\n";
         if((fabs(coef0)+fabs(coef1)+fabs(coef2))*S.delt > 1)
         {
          cout<<"************************************************FM1*****i+1,j,k="<<i<<","<<j<<","<<k<<"\n";
          cout<<"*****FM1**************************************coefs="<<coefsum<<","<<coef0<<","<<coef1<<","<<coef2<<"\n";
         }

         if(isnan(F[i+gama][j][k]) || isinf(F[i+gama][j][k]))
         {
            cout<<"***FM1*****i+1,j,k="<<i<<","<<j<<","<<k<<"\n";
            exit(0);
         }
         if(coef0>0 || coef1>0 || coef2>0) cout<<"\n Algum coef>0 \n";
         double temp=fabs(-derconv*S.delt +  U[i+gama][j][k]);
         double max_temp = 0;
               for(int idel=-2;idel<=2;idel++) 
               {
                  for(int jdel=-2;jdel<=2;jdel++)
                  {
                      for(int kdel=-2;kdel<=2;kdel++)
                      {
                        if((i+gama+idel>=1) && (j+jdel>=1) && (k+kdel>=1) && 
                        (i+gama+idel<=S.gridp[0]) && (j+jdel<=S.gridp[1]) && (k+kdel<=S.gridp[2]) )
                           if(((ISNOTAIR(flag[i+gama+idel][j+jdel][k+kdel])) || (ISNOTAIR(flag[i+gama+idel-1][j+jdel][k+kdel]))) &&
                              (fabs(U[i+gama+idel][j+jdel][k+kdel]) > max_temp))
                              max_temp = fabs(U[i+gama+idel][j+jdel][k+kdel]);
                       }
                  }
               }
               if(temp > max_temp)
               {
                   cout<<"***********************************FM1************************************ \n";
                   cout<<"**********i+1,j,k,temp,max_temp"<<i<<","<<j<<","<<k<<","<<temp<<","<<max_temp<<"\n";
               }*/
}
void Navier::CompGM1(Velocity gridCellUVW, int i, int j, int k, int gama)//calculate temporary V velocity
{
              //cout<<"\n ****************** Entrei na CompFV1 ******************************** \n";
              //cout<<"\n *****************************"<<i<<","<< j<<","<<k<<"***********************\n";
	      int alpha,beta;
	      double e0,e1,e2,c1,c2,d1,d2,f1,f2,g1,g2;
	      double coef0,coef1,coef2;
         double d2vdx2, d2vdy2, d2vdz2;
         double Gy;
         Gy = S.g[1]/S.froude;
         S.nu = 0.0;
         double h=1.0/S.gridp[0];
         DX[0] = h; DY[0]= h; DZ[0]=h;
         DX[i] =h; DY[j]=h; DZ[k]=h;
         DX[i+1] =h; DY[j+1]=h; DZ[k+1]=h;
         DX[i-1] =h; DY[j-1]=h; DZ[k-1]=h;
	      alpha = 2*gama -1;
	      beta = alpha + gama;

	      coef0 = alpha*V[i][j+gama][k];
	      if(coef0>=0)
         {
            coef0 *=(-1.0/DY[j]);
            e0= V[i][j+1-gama][k];
            coef1 = U[i+1][j][k]+U[i][j][k];
            coef2 = W[i][j][k+1]+W[i][j][k];
	      }
	     else
	     {
            coef0 *=(1.0/DY[j+alpha]);
            e0= V[i][j+gama][k];
            /*if((j<S.gridp[1]) || (alpha == -1)) 
            {
               e0 =.25*(gridCellUVW.gridV(i,j+beta,k)+gridCellUVW.gridV(i+1,j+beta,k) + 
               gridCellUVW.gridV(i,j+beta,k+1)+gridCellUVW.gridV(i+1,j+beta,k+1));
            }
            else
            {
               e0 = V[i][j+1][k];//Se o fluxo está vindo de uma ghost cell ela só pode ser inflow 
            }*/

            if(ISNOTAIR_OUTFLOW(flag[i-1][j+alpha][k]))
            {
               c1 = U[i][j+alpha][k];
            }
            else
               c1 = U[i][j][k];
            /*{
               if(j<S.gridp[1] || (alpha == -1))
               {	
                  c1 = .25*(gridCellUVW.gridU(i,j+gama,k)+gridCellUVW.gridU(i,j+beta,k) + 
                  gridCellUVW.gridU(i,j+gama,k+1)+gridCellUVW.gridU(i,j+beta,k+1));
               }
               else
               {
                  c1 = .5*(gridCellUVW.gridU(i,j+1,k)+ gridCellUVW.gridU(i,j+1,k+1));
               }
            }*/
			
            if(ISNOTAIR_OUTFLOW(flag[i+1][j+alpha][k]))
            {
               d1 = U[i+1][j+alpha][k];
            }
            else 
               d1 = U[i+1][j][k];
            /*{
               if((j<S.gridp[1]) || (alpha == -1))
               {
                  d1 =.25*(gridCellUVW.gridU(i+1,j+gama,k)+gridCellUVW.gridU(i+1,j+beta,k) + 
                  gridCellUVW.gridU(i+1,j+gama,k+1)+gridCellUVW.gridU(i+1,j+beta,k+1));
               }
               else
               {
                  d1 = .5*(gridCellUVW.gridW(i+1,j+1,k)+ gridCellUVW.gridW(i+1,j+1,k+1));
               }
            }*/
            coef1 = c1+d1;
		
            if(ISNOTAIR_OUTFLOW(flag[i][j+alpha][k-1]))
            {
               c2 = W[i][j+alpha][k];
            }
            else
               c2 = W[i][j][k];   
            /*{
               if(j<S.gridp[1] || (alpha == -1))
               {	
                  c2 = .25*(gridCellUVW.gridW(i,j+gama,k)+gridCellUVW.gridW(i,j+beta,k) + 
                  gridCellUVW.gridW(i+1,j+gama,k)+gridCellUVW.gridW(i+1,j+beta,k));
               }
               else
               {
                  c2 = .5*(gridCellUVW.gridW(i,j+1,k)+ gridCellUVW.gridW(i+1,j+1,k));
               }
            }*/
			
            if(ISNOTAIR_OUTFLOW(flag[i][j+alpha][k+1]))
            {
               d2 = W[i][j+alpha][k+1];
            }
            else 
               d2 = W[i][j][k+1];
            /*{
               if(j<S.gridp[1] || (alpha == -1))
               {   
                  d2 =.25*(gridCellUVW.gridW(i,j+gama,k+1)+gridCellUVW.gridW(i,j+beta,k+1) + 
                  gridCellUVW.gridW(i+1,j+gama,k+1)+gridCellUVW.gridW(i+1,j+beta,k+1));
               }
               else
               {
                  d2 = .5*(gridCellUVW.gridW(i,j+1,k+1)+ gridCellUVW.gridW(i+1,j+1,k+1));
               }
            }*/
            coef2 = c2 +d2;

         }
	      if((coef1)>= 0)
	      {
            coef1 *=(-1.0/(DX[i-1]+DX[i]));
            if(ISFLUID_INFLOW(flag[i-1][j][k]))
            {
               f1=e1=V[i-1][j+gama][k];
            }
            else
            {
               f1=e1=V[i][j+gama][k];
               /*f1=V[i][j+gama][k];
               e1 =.25*(gridCellUVW.gridV(i-1,j+gama,k)+gridCellUVW.gridV(i-1,j+gama,k+1) + 
               gridCellUVW.gridV(i,j+gama,k)+gridCellUVW.gridV(i,j+gama,k+1));*/
            }
            if(ISFLUID_INFLOW(flag[i+1][j][k]))
            {
               g1=V[i+1][j+gama][k];
            }
            else
            {
               g1=V[i][j+gama][k];
            }
         }
         else
         {
            coef1 *=(1.0/(DX[i+1]+DX[i]));
            if(ISFLUID_INFLOW(flag[i+1][j][k]))
            {
               g1=e1=V[i+1][j+gama][k];
            }
            else
            {	
               g1=e1=V[i][j+gama][k];
               /*g1=V[i][j+gama][k];
               if(i<S.gridp[0])
               {
                  e1 =.25*(gridCellUVW.gridV(i+1,j+gama,k)+gridCellUVW.gridV(i+1,j+gama,k+1) + 
                  gridCellUVW.gridV(i+2,j+gama,k)+gridCellUVW.gridV(i+2,j+gama,k));
               }
               else
               {
                  e1 =.5*(gridCellUVW.gridV(i+1,j+1,k)+gridCellUVW.gridV(i+1,j+1,k+1));
               }*/
            }
            if(ISFLUID_INFLOW(flag[i-1][j][k]))
            {
               f1=V[i-1][j+gama][k];
            }
            else
            {
               f1=V[i][j+gama][k];
            }
	      }
	      if((coef2)>=0)
	      {
            coef2 *=(-1.0/(DZ[k-1]+DZ[k]));
            if(ISFLUID_INFLOW(flag[i][j][k-1]))
            {
               f2=e2=V[i][j+gama][k-1];
            }
            else
            { 
              f2=e2=V[i][j+gama][k];
		        /*f2=V[i][j+gama][k];
               e2 =.25*(gridCellUVW.gridV(i,j+gama,k-1)+gridCellUVW.gridV(i+1,j+gama,k-1) + 
			     gridCellUVW.gridV(i,j+gama,k)+gridCellUVW.gridV(i+1,j+gama,k));*/
            }
            if(ISFLUID_INFLOW(flag[i][j][k+1]))
            {
               g2=V[i][j+gama][k+1];
            }
            else
            {
               g2=V[i][j+gama][k];
            }
         }
         else
         {
            coef2 *=(1.0/(DZ[k+1]+DZ[k]));
            if(ISFLUID_INFLOW(flag[i][j][k+1]))
            {
               g2=e2=V[i][j+gama][k+1];
            }
            else
            {	
               g2=e2=V[i][j+gama][k];
               /*g2=V[i][j+gama][k];
               if(k<S.gridp[2])
               {
                  e2 =.25*(gridCellUVW.gridV(i,j+gama,k+1)+gridCellUVW.gridV(i+1,j+gama,k+1) + 
			        gridCellUVW.gridV(i,j+gama,k+2)+gridCellUVW.gridV(i+1,j+gama,k+2));
               }
               else
               {
                  e2 =.5*(gridCellUVW.gridV(i,j+1,k+1)+gridCellUVW.gridV(i+1,j+1,k+1));
               }*/
            }
            if(ISFLUID_INFLOW(flag[i][j][k+1]))
            {
               f2=V[i][j+gama][k-1];
            }
            else
            {
		         f2=V[i][j+gama][k];
            }
         }
         double derconv= coef0 *(e0-V[i][j+gama][k])+ coef1*(e1-V[i][j+gama][k]) + coef2*(e2-V[i][j+gama][k]);

         d2vdy2=DD(V[i][j][k],V[i][j+gama][k],V[i][j+1][k],1,j+gama);
         d2vdx2=DDS(f1,V[i][j+gama][k],g1,0,i);
         d2vdz2=DDS(f2,V[i][j+gama][k],g2,2,k);

         G[i][j+gama][k] = S.nu*(d2vdx2+d2vdy2+d2vdz2) - derconv + Gy;
         /*
         double coef_sum =(fabs(coef0)+fabs(coef1)+fabs(coef2))*S.delt;
         if(fabs(coef0) > (UVWmax/h)) cout<<"coef0 > UVWMax="<<coef0<<","<<UVWmax/h<<"\n";
         if(fabs(coef1) > (UVWmax/h)) cout<<"coef1 > UVWMax="<<coef1<<","<<UVWmax/h<<"\n";
         if(fabs(coef2) > (UVWmax/h)) cout<<"coef2 > UVWMax="<<coef2<<","<<UVWmax/h<<"\n";
         if((fabs(coef0)+fabs(coef1)+fabs(coef2))*S.delt > 1)
         {
          cout<<"************************************************GM1*****i+1,j,k="<<i<<","<<j<<","<<k<<"\n";
          cout<<"*****GM1**************************************coefs="<<coef_sum<<","<<coef0<<","<<coef1<<","<<coef2<<"\n";
         }

         if(isnan(G[i][j+gama][k]) || isinf(G[i][j+gama][k]))
         {
            cout<<"*****GM1*****i+1,j,k="<<i<<","<<j<<","<<k<<"\n";
            exit(0);
         }
         if(coef0>0 || coef1>0 || coef2>0) cout<<"\n Algum coef>0 \n";
         double temp=fabs(-derconv*S.delt +  V[i][j+gama][k]);
         double max_temp = 0;
         if(i>1 && i<S.gridp[0] &&  j>1 && j<S.gridp[1] && k>1 && k<S.gridp[2])
         {
            for(int idel=-2;idel<=2;idel++) 
            {
               for(int jdel=-2;jdel<=2;jdel++)
               {
                  for(int kdel=-2;kdel<=2;kdel++)
                  {
                     if((i+idel>=1) && (j+gama+jdel>=1) && (k+kdel>=1) && 
                     (i+idel<=S.gridp[0]) && (j+gama+jdel<=S.gridp[1]) && (k+kdel<=S.gridp[2]) )
                        if(((ISNOTAIR(flag[i+idel][j+gama+jdel][k+kdel])) || (ISNOTAIR(flag[i+idel][j+gama+jdel-1][k+kdel]))) &&
                          (fabs(V[i+idel][j+gama+jdel][k+kdel]) > max_temp))
                           max_temp = fabs(V[i+idel][j+gama+jdel][k+kdel]);
                  }
               }
            }
            if(temp > max_temp)
            {
                cout<<"*********************************GM1************************************** \n";
                cout<<"**********i+1,j,k,temp,max_temp"<<i<<","<<j<<","<<k<<","<<temp<<","<<max_temp<<"\n";
            }
         }*/
}

void Navier::CompHM1(Velocity gridCellUVW, int i, int j, int k, int gama)//calculate temporary W velocity
{
   //cout<<"\n ****************** Entrei na CompFV1 ******************************** \n";
   //cout<<"\n *****************************"<<i<<","<< j<<","<<k<<"***********************\n";
   int alpha,beta;
   double e0,e1,e2,c1,c2,d1,d2,f1,f2,g1,g2;
   double coef0,coef1,coef2;
   double d2wdx2, d2wdy2, d2wdz2;
   double Gz;
   Gz = S.g[2]/S.froude;
   S.nu = 0.0;
   double h=1.0/S.gridp[0];
   DX[0] =h; DY[0]=h; DZ[0]=h;
   DX[i] =h; DY[j]=h; DZ[k]=h;
   DX[i+1] =h; DY[j+1]=h; DZ[k+1]=h;
   DX[i-1] =h; DY[j-1]=h; DZ[k-1]=h;
   
   alpha = 2*gama -1;
   beta = alpha + gama;

   coef0 = alpha*W[i][j][k+gama];
   if(coef0>=0)
   {
      coef0 *=(-1.0/DZ[j]);
      e0= W[i][j][k+1-gama];
      coef1 = U[i+1][j][k]+U[i][j][k];
      coef2 = V[i][j+1][k]+V[i][j][k];
   }
   else
   {
      coef0 *=(1.0/DZ[j+alpha]);
      e0= W[i][j][k+gama];
      /*if((k<S.gridp[2]) || (alpha == -1)) 
      {
         e0 =.25*(gridCellUVW.gridW(i,j,k+beta)+gridCellUVW.gridW(i+1,j,k+beta) + 
         gridCellUVW.gridW(i,j+1,k+beta)+gridCellUVW.gridW(i+1,j+1,k+beta));
      }
      else
      {
         e0 = W[i][j][k+1];//Se o fluxo está vindo de uma ghost cell ela só pode ser inflow 
      }*/

      if(ISNOTAIR_OUTFLOW(flag[i-1][j][k+alpha]))
      {
         c1 = U[i][j][k+alpha];
      }
      else 
         c1 = U[i][j][k];
      /*{
         if(k<S.gridp[2] || (alpha == -1))
         {	
            c1 = .25*(gridCellUVW.gridU(i,j,k+gama)+gridCellUVW.gridU(i,j,k+beta) + 
			   gridCellUVW.gridU(i,j+1,k+gama)+gridCellUVW.gridU(i,j+1,k+beta));
         }
         else
         {
            c1 = .5*(gridCellUVW.gridU(i,j,k+1)+ gridCellUVW.gridU(i,j+1,k+1));
         }
      }*/
			
      if(ISNOTAIR_OUTFLOW(flag[i+1][j][k+alpha]))
      {
         d1 = U[i+1][j][k+alpha];
      }
      else 
         d1 = U[i+1][j][k];
      /*{
         if(k<S.gridp[2] || (alpha == -1))
         {
            d1 =.25*(gridCellUVW.gridU(i+1,j,k+gama)+gridCellUVW.gridU(i+1,j,k+beta) + 
			     gridCellUVW.gridU(i+1,j+1,k+gama)+gridCellUVW.gridU(i+1,j+1,k+beta));
         }
         else
         {
			   d1 = .5*(gridCellUVW.gridV(i+1,j,k+1)+ gridCellUVW.gridV(i+1,j+1,k+1));
         }
      }*/
      coef1 = c1+d1;
		
      if(ISNOTAIR_OUTFLOW(flag[i][j-1][k+alpha]))
      {
         c2 = V[i][j][k+alpha];
      }
      else 
         c2 = V[i][j][k];
      /*{
         if(k<S.gridp[2] || (alpha == -1))
         {	
            c2 = .25*(gridCellUVW.gridV(i,j,k+gama)+gridCellUVW.gridV(i,j,k+beta) + 
			   gridCellUVW.gridV(i+1,j,k+gama)+gridCellUVW.gridV(i+1,j,k+beta));
         }
         else
         {
            c2 = .5*(gridCellUVW.gridV(i,j,k+1)+ gridCellUVW.gridV(i+1,j,k+1));
         }
      }*/
			
      if(ISNOTAIR_OUTFLOW(flag[i][j+1][k+alpha]))
      {
         d2 = V[i][j+1][k+alpha];
      }
      else 
         d2 = V[i][j+1][k];
      /*{
         if(k<S.gridp[2] || (alpha == -1))
         {
            d2 =.25*(gridCellUVW.gridV(i,j+1,k+gama)+gridCellUVW.gridV(i,j+1,k+beta) + 
			     gridCellUVW.gridV(i+1,j+1,k+gama)+gridCellUVW.gridV(i+1,j+1,k+beta));
         }
         else
         {
            d2 = .5*(gridCellUVW.gridV(i,j+1,k+1)+ gridCellUVW.gridV(i+1,j+1,k+1));
         }
      }*/
      coef2 = c2 +d2;
   }
   if((coef1)>= 0)
   {
      coef1 *=(-1.0/(DX[i-1]+DX[i]));
      if(ISFLUID_INFLOW(flag[i-1][j][k]))
      {
         f1=e1=W[i-1][j][k+gama];
      }
      else
      {
         f1=e1=W[i][j][k+gama];
         /*f1=W[i][j][k+gama];
         e1 =.25*(gridCellUVW.gridW(i-1,j,k+gama)+gridCellUVW.gridW(i-1,j+1,k+gama) + 
         gridCellUVW.gridW(i,j,k+gama)+gridCellUVW.gridW(i,j+1,k+gama));*/
      }
      if(ISFLUID_INFLOW(flag[i+1][j][k]))
      {
         g1=W[i+1][j][k+gama];
      }
      else
      {
         g1=W[i][j][k+gama];
      }
   }
   else
   {
      coef1 *=(1.0/(DX[i+1]+DX[i]));
      if(ISFLUID_INFLOW(flag[i+1][j][k]))
      {
         g1=e1=W[i+1][j][k+gama];
      }
      else
      {	
          g1=e1=W[i][j][k+gama];
         /*g1=W[i][j][k+gama];
         if(i<S.gridp[0])
         {
            e1 =.25*(gridCellUVW.gridW(i+1,j,k+gama)+gridCellUVW.gridW(i+1,j+1,k+gama) + 
            gridCellUVW.gridW(i+2,j,k+gama)+gridCellUVW.gridW(i+2,j,k+gama));
         }
         else
         {
            e1 =.5*(gridCellUVW.gridW(i+1,j,k+1)+gridCellUVW.gridW(i+1,j,k+1+1));
         }*/
      }
      if(ISFLUID_INFLOW(flag[i-1][j][k]))
      {
         f1=W[i-1][j][k+gama];
      }
      else
      {
         f1=W[i][j][k+gama];
         }}

   if((coef2)>=0)
   {
      coef2 *=(-1.0/(DY[j-1]+DY[j]));
      if(ISFLUID_INFLOW(flag[i][j-1][k]))
         {
            f2=e2=W[i][j-1][k+gama];
         }
      else
      {
         f2=e2=W[i][j][k+gama];
         /*f2=W[i][j][k+gama];
         e2 =.25*(gridCellUVW.gridW(i,j-1,k+gama)+gridCellUVW.gridW(i+1,j-1,k+gama) + 
         gridCellUVW.gridW(i,j,k+gama)+gridCellUVW.gridW(i+1,j,k+gama));
         */
      }
      if(ISFLUID_INFLOW(flag[i][j+1][k]))
      {
         g2=W[i][j+1][k+gama];
      }
      else
      {
        g2=W[i][j][k+gama];
      }
   }
   else
{
      coef2 *=(1.0/(DY[j+1]+DY[j]));
      if(ISFLUID_INFLOW(flag[i][j+1][k]))
      {
         g2=e2=W[i][j+1][k+gama];
      }
      else
      {	
         g2=e2=W[i][j][k+gama];
         /*g2=W[i][j][k+gama];
         if(j<S.gridp[1])
         {
            e2 =.25*(gridCellUVW.gridW(i,j+1,k+gama)+gridCellUVW.gridW(i+1,j+1,k+gama) + 
           gridCellUVW.gridW(i,j+2,k+gama)+gridCellUVW.gridW(i+1,j+2,k+gama));
         }
         else
         {
         e2 =.5*(gridCellUVW.gridW(i,j+1,k+1)+gridCellUVW.gridW(i+1,j+1,k+1));
         }*/
      }
      if(ISFLUID_INFLOW(flag[i][j+1][k]))
      {
         f2=W[i][j-1][k+gama];
      }
      else
      {
        f2=W[i][j][k+gama];
      }
   }
   double derconv= coef0 *(e0-W[i][j][k+gama])+ coef1*(e1-W[i][j][k+gama]) + coef2*(e2-W[i][j][k+gama]);

   d2wdz2=DD(W[i][j][k],W[i][j][k+gama],W[i][j][k+1],2,k+gama);
   d2wdx2=DDS(f1,W[i][j][k+gama],g1,0,i);
   d2wdy2=DDS(f2,W[i][j][k+gama],g2,1,j);


   H[i][j][k+gama] = S.nu*(d2wdx2+d2wdy2+d2wdz2) - derconv + Gz; 
   
   /*if(fabs(coef0) > (UVWmax/h)) cout<<"coef0 > UVWMax="<<coef0<<","<<UVWmax/h<<"\n";
   if(fabs(coef1) > (UVWmax/h)) cout<<"coef1 > UVWMax="<<coef1<<","<<UVWmax/h<<"\n";
   if(fabs(coef2) > (UVWmax/h)) cout<<"coef2 > UVWMax="<<coef2<<","<<UVWmax/h<<"\n";
   double coefsum =(fabs(coef0)+fabs(coef1)+fabs(coef2))*S.delt;
   if((fabs(coef0)+fabs(coef1)+fabs(coef2))*S.delt > 1)
   {
    cout<<"************************************************HM1*****i+1,j,k="<<i<<","<<j<<","<<k<<"\n";
    cout<<"*****HM1**************************************coefs="<<coefsum<<","<<coef0<<","<<coef1<<","<<coef2<<"\n";
   }
   if(isnan(H[i][j][k+gama]) || isinf(H[i][j][k+gama]))
               {
                  cout<<"*****HM1*****i+1,j,k="<<i<<","<<j<<","<<k<<"\n";
                  exit(0);
               }
   if(coef0>0 || coef1>0 || coef2>0) cout<<"\n Algum coef>0 \n";
   double temp=fabs(-derconv*S.delt +  W[i][j][k+gama]);
   double max_temp = 0;
   if(i>1 && i<S.gridp[0] &&  j>1 && j<S.gridp[1] && k>1 && k<S.gridp[2])
   {
      for(int idel=-2;idel<=2;idel++) 
      {
         for(int jdel=-2;jdel<=2;jdel++)
         {
             for(int kdel=-2;kdel<=2;kdel++)
             {
                if((i+idel>=1) && (j+jdel>=1) && (k+gama+kdel>=1) && 
                   (i+idel<=S.gridp[0]) && (j+jdel<=S.gridp[1]) && (k+gama+kdel<=S.gridp[2]) )
                if(((ISNOTAIR(flag[i+idel][j+jdel][k+gama+kdel])) || (ISNOTAIR(flag[i+idel][j+jdel][k+gama+kdel-1]))) &&
                    (fabs(W[i+idel][j+jdel][k+gama+kdel]) > max_temp))
                       max_temp = fabs(W[i+idel][j+jdel][k+gama+kdel]);
              }
         }
      }
      if(temp > max_temp)
      {
          cout<<"********************************HM1*************************************** \n";
          cout<<"**********i+1,j,k,temp,max_temp"<<i<<","<<j<<","<<k<<","<<temp<<","<<max_temp<<"\n";
      }
   }*/
}
//Compute auxiliary-velocities without free-boundaries ##
void Navier::CompFGH(Velocity gridCellUVW)//calculate temporary velocity fields (in space) due to chorin projection method (see User's Guide) 
{
   int i,j,k;
   double Fmax = 0, Gmax = 0, Hmax = 0;
   
   int d1[6]={-1,1,0,0,0,0};
   int d2[6]={0,0,-1,1,0,0};
   int d3[6]={0,0,0,0,-1,1};
   int factor0[6]={1,1,1,1,1,1};
   int factor1[6]={0,0,0,0,0,0};
   int mn[6]={1,0,3,2,5,4};
   //<=5 point stencils (VONOS, SMART, HLPA, QUICK, DC) 
   JLOOP 
      KLOOP 
         for (i=1; i<=S.gridp[0]; i++) 
         {
            IFFLUID(flag[i][j][k])
            {  
               if(flag[i][j][k]%OFFSET_LIMIT != 0)//tratando a c�lulas slip's que possuem mais de uma c�lula livre (air ou fluid)
               {
                  for(int n=0;n<=5;n++)
                  {
                     if(flag[i+d1[n]][j+d2[n]][k+d3[n]]== SLIP_ZERO_TWO_FREE_FACES)
                     {
                        int m=mn[n];
                        factor0[n] = 0;
                        //factor1[n] = 0;
                        factor0[m] = 0;
                        factor1[m] = 0;
                        Update_Slip_Outflow_Cells(i,j,k,i+d1[n],j+d2[n],k+d3[n],factor0,factor1);
                        factor0[m] = 1;
                        factor1[m] = 0;
                        factor0[n] = 1;
                        //factor1[n] = 1;
                     }
                     else
                     {
                        if(flag[i+d1[n]][j+d2[n]][k+d3[n]]== OUTFLOW_ZERO_TWO_FREE_FACES)
                        {
                           //int m=mn[n];
                           factor0[n] = 0;
                           factor1[n] = 1;
                           Update_Slip_Outflow_Cells(i,j,k,i+d1[n],j+d2[n],k+d3[n],factor0, factor1);
                           factor0[n] = 1;
                           factor1[n] = 0;
                        }
                     }
                  }
               }
               if(ISNOTGIVEN(flag[i-1][j][k]))
               {
                  if(flag[i-1][j][k] > AIR_FLUID_END) 
                      CompF(i,j,k);
                  else
                    CompFM1(gridCellUVW,i,j,k,0);
               }
               if(flag[i+1][j][k]  <= AIR_FLUID_END)
                  CompFM1(gridCellUVW,i,j,k,1);
               if(fabs(F[i][j][k]) > fabs(Fmax)) Fmax = F[i][j][k];
               
               if(ISNOTGIVEN(flag[i][j-1][k])) 
               {
                  if(flag[i][j-1][k] > AIR_FLUID_END)                
                     CompG(i,j,k); 
                 else
                  CompGM1(gridCellUVW,i,j,k,0);
               }
               if(flag[i][j+1][k] <= AIR_FLUID_END)
                  CompGM1(gridCellUVW, i,j,k,1);
               if(fabs(G[i][j][k]) > fabs(Gmax)) Gmax = G[i][j][k];
               
               if(ISNOTGIVEN(flag[i][j][k-1])) 
               {
                  if(flag[i][j][k-1] > AIR_FLUID_END)  
                     CompH(i,j,k);
                  else
                  CompHM1(gridCellUVW,i,j,k,0);
               }
               if(flag[i][j][k+1] <= AIR_FLUID_END)
                  CompHM1(gridCellUVW, i,j,k,1);
               if(fabs(H[i][j][k]) > fabs(Hmax)) Hmax = H[i][j][k];
                  
            }
         }
}

//time handling for auxiliary-velocities
void Navier::CompTUVWfromFGH(int timestepmethod) 
{  
   int i,j,k;
   double uvwmax = 0;
   int imax, jmax, kmax;
    double uvwmaxantes = 0;
   int imaxantes, jmaxantes, kmaxantes;
   //cout<<"/n Entrando na TUVW\n";
   	//##############################
      //##     EULER 1st ORDER      ##
      //##############################
      JLOOP 
         KLOOP 
            ILOOP 
               //if(flag[i][j-1][k]!=INFLOW)   
               IFFLUID(flag[i][j][k])
               {
                   if(fabs(U[i][j][k])>fabs(uvwmaxantes))
                  {
                     uvwmaxantes = U[i][j][k];
                     imaxantes = i; jmaxantes = j; kmaxantes = k;
                  }
                   /*if((i<=3) && (j>=10) && (j<=14) && (k>=5) && (k<=9))
                     {
                        cout<<"***************************************************************************\n";
                        cout<<"/n i,j,k,U antes ="<<i<<","<<j<<","<<k<<","<<U[i][j][k]<<"\n";
                     }*/
                  if(ISNOTGIVEN(flag[i-1][j][k]))
                     U[i][j][k]+=S.delt*F[i][j][k];
                  if(flag[i+1][j][k] <= AIR_FLUID_END)
                     U[i+1][j][k]+=S.delt*F[i+1][j][k];
                    /* if((i<=3) && (j>=10) && (j<=14) && (k>=5) && (k<=9))
                     {
                        cout<<"***************************************************************************\n";
                        cout<<"/n i,j,k,U ="<<i<<","<<j<<","<<k<<","<<U[i][j][k]<<"\n";
                        cout<<"/n i,j,k,F="<<i<<","<<j<<","<<k<<","<<F[i][j][k]<<"\n";
                        cout<<"/n i,j,k,F+1 ="<<i<<","<<j<<","<<k<<","<<F[i+1][j][k]<<"\n";
                        
                     }*/
                 
                  if(U[i][j][k]>fabs(uvwmax))
                  {
                     uvwmax = U[i][j][k];
                     imax = i; jmax = j; kmax = k;
                  }
               }
      ILOOP 
         KLOOP 
            JLOOP            
               //if(flag[i][j-1][k]!=INFLOW)
               IFFLUID(flag[i][j][k])
               {
                   if(fabs(V[i][j][k])>fabs(uvwmaxantes))
                  {
                     uvwmaxantes = V[i][j][k];
                     imaxantes = i; jmaxantes = j; kmaxantes = k;
                  }
                  if(ISNOTGIVEN(flag[i][j-1][k]))
                     V[i][j][k]+=S.delt*G[i][j][k];
                  if(flag[i][j+1][k] <= AIR_FLUID_END)
                     V[i][j+1][k]+=S.delt*G[i][j+1][k];                  
                  /*if((i<=3) && (j>=10) && (j<=14) && (k>=5) && (k<=9))
                     {
                        cout<<"***************************************************************************\n";
                        cout<<"/n i,j,k,V ="<<i<<","<<j<<","<<k<<","<<V[i][j][k]<<"\n";
                        cout<<"/n i,j,k,G ="<<i<<","<<j<<","<<k<<","<<G[i][j][k]<<"\n";
                        cout<<"/n i,j,k,G+1 ="<<i<<","<<j<<","<<k<<","<<G[i][j+1][k]<<"\n";
                        
                     }*/

                  if(V[i][j][k]>fabs(uvwmax))
                  {
                     uvwmax = V[i][j][k];
                     imax = i; jmax = j; kmax = k;
                  }
               }
      ILOOP 
         JLOOP 
            KLOOP
               //if(flag[i][j-1][k]!=INFLOW)
               IFFLUID(flag[i][j][k])
               {
                   if(fabs(W[i][j][k])>fabs(uvwmaxantes))
                  {
                     uvwmaxantes = W[i][j][k];
                     imaxantes = i; jmaxantes = j; kmaxantes = k;
                  }
                  if(ISNOTGIVEN(flag[i][j][k-1]))
                     W[i][j][k]+=S.delt*H[i][j][k];
                  if(flag[i][j][k+1] <= AIR_FLUID_END)
                     W[i][j][k+1]+=S.delt*H[i][j][k+1];                  
                  /*if((i<=3) && (j>=9) && (j<=13) && (k<=3))
                     {
                        cout<<"***************************************************************************\n";
                        cout<<"/n i,j,k,W ="<<i<<","<<j<<","<<k<<","<<W[i][j][k]<<"\n";
                        cout<<"/n i,j,k,H="<<i<<","<<j<<","<<k<<","<<H[i][j][k]<<"\n";
                        cout<<"/n i,j,k,H+1 ="<<i<<","<<j<<","<<k<<","<<H[i][j][k+1]<<"\n";
                        
                     }*/

                  if(fabs(W[i][j][k])>fabs(uvwmax))
                  {
                     uvwmax = W[i][j][k];
                     imax = i; jmax = j; kmax = k;
                  }
               }   
}

//Compute right hand side for poisson equation
void Navier::CompRHS() 
{
   int i, j, k;
   IJKLOOP    
      //IFFLUID(flag[i][j][k])
      if(ISFLUID_OUTFLOW(flag[i][j][k]))
         RHS[i][j][k]=  ((U[i+1][j][k]-U[i][j][k])/(DX[i]) +
                        (V[i][j+1][k]-V[i][j][k])/(DY[j])+
                        (W[i][j][k+1]-W[i][j][k])/(DZ[k])) /*/S.delt */ ;
      
}

//calc new velocities by pressure-gradient correction  
//c�lcula a velocidade final, ap�s o c�lculo da press�o que satisfaz a condi��o de diverg�ncia livre
void Navier::AdapUVW() 
{
   int i,j,k;
   /*double uvmaxproj = 0.0;
   int imax,jmax,kmax,control;
      IJKLOOP
   {
      //double div = 0.0;
      IFFLUID(flag[i][j][k])
      {
         if((j>=1)&&(j<=3) && (k>=6) && (k<=8) && (i>=10) && (i<=12) && (num_iter_Navier > 40))
         {
            cout<<"i,j="<<i<<","<<j<<","<<k<<","<<U[i][j][k]<<","<<V[i][j][k]<<","<<W[i][j][k]<<"\n";
            cout<<"i,j="<<i<<","<<j<<","<<U[i+1][j][k]<<","<<V[i][j+1][k]<<","<<W[i][j][k+1]<<"\n";
            //cout<<"i,j="<<i<<","<<j<<","<<P[i][j][k]<<","<<P[i][j+1][k]<<"\n\n";
         }
         if(fabs(U[i][j][k]) > fabs(uvmaxproj)) {uvmaxproj = U[i][j][k]; imax =i;jmax=j; kmax=k; control = 1;}
         if(fabs(V[i][j][k]) > fabs(uvmaxproj)) {uvmaxproj = V[i][j][k];  imax =i;jmax=j; kmax=k; control =2;}
         if(fabs(W[i][j][k]) > fabs(uvmaxproj)) {uvmaxproj = W[i][j][k];  imax =i;jmax=j; kmax=k; control =3;}
         IFFLUID(flag[i+1][j][k])
            IFFLUID(flag[i][j+1][k])
               IFFLUID(flag[i][j][k+1])
               {   
                  div = U[i+1][j][k] + V[i][j+1][k] + W[i][j][k+1] - (U[i][j][k] + V[i][j][k] + W[i][j][k]);
                  if (fabs(div) > fabs(divMax))
                     divMax = div;                  
               }
      }
   }
   cout<<"\n**********************Maximal Value Antes="<<uvmaxproj<<"\n";
   cout<<"i,j,k,control ="<<imax<<","<<jmax<<","<<kmax<<","<<control<<"\n";*/

   JLOOP 
      KLOOP 
         ILOOP 
            IFFLUID(flag[i][j][k])             
            {
               if(ISFLUID_OUTFLOW(flag[i-1][j][k]))
                  U[i][j][k]-=/*S.delt* */ (P[i][j][k]-P[i-1][j][k])/DXM[i+1];//j� foi multiplicada por S.delt antes da chamada da fun��o   
               else
                  if(ISAIR(flag[i-1][j][k]))
                     U[i][j][k]-=/*S.delt* */ P[i][j][k]/DXM[i+1];
               
               if(ISAIR(flag[i+1][j][k]))
                  U[i+1][j][k]+=/*S.delt* */ P[i][j][k]/DXM[i+1];
   
               if(ISFLUID_OUTFLOW(flag[i][j-1][k]))
                  V[i][j][k]-=/*S.delt* */ (P[i][j][k]-P[i][j-1][k])/DYM[j+1];
               else
                  if(ISAIR(flag[i][j-1][k]))
                     V[i][j][k]-=/*S.delt* */ (P[i][j][k])/DYM[j+1];
               
               if(ISAIR(flag[i][j+1][k]))
                  V[i][j+1][k]+=/*S.delt* */ P[i][j][k]/DYM[j+1];
  
               if(ISFLUID_OUTFLOW(flag[i][j][k-1]))
                  W[i][j][k]-=/*S.delt* */ (P[i][j][k]-P[i][j][k-1])/DZM[k+1];
               else
                  if(ISAIR(flag[i][j][k-1]))
                     W[i][j][k]-=/*S.delt* */ P[i][j][k]/DZM[k+1];
               
               if(ISAIR(flag[i][j][k+1]))
                  W[i][j][k+1]+=/*S.delt* */ P[i][j][k]/DZM[k+1];
            }
            
  // double divMax = 0.0;
  // uvmaxproj = 0.0;
   
  /* IJKLOOP
   {
      //double div = 0.0;
      IFFLUID(flag[i][j][k])
      {
        if(j==1)
         {
            cout<<"i,j="<<i<<","<<j<<","<<k<<","<<U[i][j][k]<<","<<V[i][j][k]<<","<<W[i][j][k]<<"\n";
            cout<<"i,j="<<i<<","<<j<<","<<U[i+1][j][k]<<","<<V[i][j+1][k]<<","<<W[i][j][k+1]<<"\n";
            cout<<"i,j="<<i<<","<<j<<","<<P[i][j][k]<<","<<P[i][j+1][k]<<"\n\n";
         }
         //if(fabs(U[i][j][k]) > fabs(uvmaxproj)) {uvmaxproj = U[i][j][k]; imax =i;jmax=j; kmax=k; control = 1;}
         //if(fabs(V[i][j][k]) > fabs(uvmaxproj)) {uvmaxproj = V[i][j][k];  imax =i;jmax=j; kmax=k; control =2;}
         //if(fabs(W[i][j][k]) > fabs(uvmaxproj)) {uvmaxproj = W[i][j][k];  imax =i;jmax=j; kmax=k; control =3;}
         IFFLUID(flag[i+1][j][k])
            IFFLUID(flag[i][j+1][k])
               IFFLUID(flag[i][j][k+1])
               {   
                  div = U[i+1][j][k] + V[i][j+1][k] + W[i][j][k+1] - (U[i][j][k] + V[i][j][k] + W[i][j][k]);
                  if (fabs(div) > fabs(divMax))
                     divMax = div;                  
               }
      }
   } */
   //cout<<"\n**********************Maximal Value ="<<uvmaxproj<<"\n";
   //cout<<"i,j,k,control ="<<imax<<","<<jmax<<","<<kmax<<","<<control<<"\n";
  //cout<<"/n DIVMAX---------------->"<<divMax<<"/n";*/
}

//Time-Step Control
double Navier::TimeStep()
{
   double deltneu=S.deltmax ; //S.deltmax � definido no arquivo .nav neste caso = 0.01
   double  xmin=1e15,ymin=1e15,zmin=1e15 ;
   double uvwmax = 0.0;
   int i,j,k;
   if(xmin==1e15 && ymin==1e15 && zmin==1e15)
   {                   
      IJKLOOP 
      {
         xmin=min(xmin,DX[i]); 
         ymin=min(ymin,DY[j]); 
         zmin=min(zmin,DZ[k]);
         if(fabs(U[i][j][k]) > uvwmax)  uvwmax = fabs(U[i][j][k]);
         if(fabs(V[i][j][k]) > uvwmax)  uvwmax = fabs(V[i][j][k]);
         if(fabs(W[i][j][k]) > uvwmax)  uvwmax = fabs(W[i][j][k]);
      } 
   }
   double xyzmin = min(min(xmin,ymin),zmin);
   // local CFL-condition including volume forces
   IJKLOOP
   {  //S.tfconv � definido no arquivo typen.cpp 25 neste caso = 0.1; S.froude = 1.0
      deltneu=min(deltneu, S.tfconv*2/(fabs(U[i][j][k])/DXM[i]+
				       sqrt( (U[i][j][k]*U[i][j][k])/(DXM[i]*DXM[i])+4*fabs(S.g[0])/(S.froude*DXM[i])) +1E-15));
      deltneu=min(deltneu, S.tfconv*2/(fabs(V[i][j][k])/DYM[j]+
				       sqrt( (V[i][j][k]*V[i][j][k])/(DYM[j]*DYM[j])+4*fabs(S.g[1])/(S.froude*DYM[j])) +1E-15));
      deltneu=min(deltneu, S.tfconv*2/(fabs(W[i][j][k])/DZM[k]+
				       sqrt( (W[i][j][k]*W[i][j][k])/(DZM[k]*DZM[k])+4*fabs(S.g[2])/(S.froude*DZM[k])) +1E-15));
    }  
      double lim_temp = xyzmin/(3.0*uvwmax);
      deltneu=min(deltneu,lim_temp);
      UVWmax = uvwmax;
      cout<<"******************************************UVWmax ="<<UVWmax<<"\n";
    // Viscosity 
    //tfdiff = 0.3;
    deltneu=min(deltneu,S.tfdiff*S.re*0.5/(1./(xmin*xmin)+1./(ymin*ymin)+1./(zmin*zmin))); 
    cout<<"\n deltneu,xyzmin ="<<deltneu<<","<<xyzmin<<"\n";
   if(deltneu<1E-15) 
   {
      printf("Maybe ERROR in parameters: Time-Step = 0.0 ! \n");
      exit(1);    
   } 
   //deltneu = deltneu/2;
   return deltneu;
}

Navier::~Navier() { }
//define a velocidade nas faces "posteriore" em aberto, das c�lulas da borda, essa defini��o � feita para facilitar
//a defini��o inicial da velocidade nas novas c�lulas Fluid (AirToFluid)
void Navier::Complete_Cell(int i,int j, int k, int freefacespar)
{
   double  tempMx=0, tempMy=0, tempMz=0;
   int cardx=0, cardy=0, cardz=0;
   double deltax,deltay,deltaz;
   double excessxby2,excessyby2,excesszby2,excessby3;
   double mdeltax, mdeltay,mdeltaz;

   switch(freefacespar)
   {
      case 0:
      {
         if(flag[i-1][j][k]>= AIR_FLUID_END)//se True extrapola o valor de Ui+1 a partir Ui-1 e Ui
         {
            tempMx += 2* U[i][j][k] - U[i-1][j][k];
            cardx ++;
            if((flag[i-1][j+1][k]>= AIR_FLUID_END)|| (flag[i-1][j][k]>= AIR_FLUID_END) )
            {
               tempMy += V[i-1][j+1][k];
               cardy ++;
            }
            if((flag[i-1][j][k+1]>= AIR_FLUID_END)|| (flag[i-1][j][k]>= AIR_FLUID_END))
            {
               tempMz += W[i-1][j][k+1];
               cardz ++;
            }
         }
         else//se i-1 n�o estiver definido, ao inv�s de extrapolar repito o valor de Ui em Ui+1
         {
            tempMx += U[i][j][k];
            cardx ++;
         }

         if(flag[i][j-1][k]>= AIR_FLUID_END)//verificando a outra face vizinha anterior 
         {
            tempMy += 2* V[i][j][k] - V[i][j-1][k];
            cardy ++;
            if((flag[i+1][j-1][k]>= AIR_FLUID_END)|| (flag[i][j-1][k]>= AIR_FLUID_END))
            {
               tempMx += U[i+1][j-1][k];
               cardx ++;
            }
            if((flag[i][j-1][k+1]>= AIR_FLUID_END)|| (flag[i][j-1][k]>= AIR_FLUID_END))
            {
               tempMz += W[i][j-1][k+1];
               cardz ++;
            }
         }
         else
         {
            tempMy += V[i][j][k];
            cardy ++;
         }
         if(flag[i][j][k-1]>= AIR_FLUID_END)
         {
            tempMz += 2* W[i][j][k] - W[i][j][k-1];
            cardz ++;
            if((flag[i+1][j][k-1]>= AIR_FLUID_END)||(flag[i][j][k-1]>= AIR_FLUID_END))
            {
               tempMx += U[i+1][j][k-1];
               cardx++;
            }
            if((flag[i][j+1][k-1]>= AIR_FLUID_END)||(flag[i][j][k-1]>= AIR_FLUID_END))
            {
               tempMy += V[i][j+1][k-1];
               cardy++;
            }
         }
         else
         {
            tempMz += W[i][j][k];
            cardz ++;
         }
         tempMx/=cardx;
         tempMy/=cardy;
         tempMz/=cardz;
         deltax = tempMx - U[i][j][k];//aproxima��o de du/dx na c�lula
         deltay = tempMy - V[i][j][k];//aproxima��o de dv/dy na c�lula
         deltaz = tempMz - W[i][j][k];//aproxima��o de dw/dz na c�lula
         excessby3 = (- deltaz - deltax - deltay)/3;//calculo do excesso (divergencia deve ser igual a zero)
         U[i+1][j][k] = tempMx + excessby3;//distribuindo igualmente o excesso nas 3 dire��es
         V[i][j+1][k] = tempMy + excessby3;//   "
         W[i][j][k+1] = tempMz + excessby3;//   "   
         //cout<<"CompCell 1";
         return;
         }      
      case 1:
      {
         if((flag[i+1][j][k+1]>= AIR_FLUID_END)||(flag[i][j][k+1]>= AIR_FLUID_END))
         {
            tempMx += U[i+1][j][k+1];
            cardx++;
         }
         if((flag[i][j+1][k+1]>= AIR_FLUID_END)||(flag[i][j][k+1]>= AIR_FLUID_END))
         {
            tempMy += V[i][j+1][k+1];
            cardy++;
         }
         if(flag[i-1][j][k]>= AIR_FLUID_END)
         {
            tempMx += 2* U[i][j][k] - U[i-1][j][k];
            cardx++;
            if((flag[i-1][j+1][k]>= AIR_FLUID_END)||(flag[i-1][j][k]>= AIR_FLUID_END) )
            {
               tempMy += V[i-1][j+1][k];
               cardy++;
            }
         }
         else
         {
            tempMx += U[i][j][k];
            cardx ++;
         }
         if(flag[i][j-1][k]>= AIR_FLUID_END)
         {
            tempMy += 2* V[i][j][k] - V[i][j-1][k];
            cardy++;
            if((flag[i+1][j-1][k]>= AIR_FLUID_END)|| (flag[i][j-1][k]>= AIR_FLUID_END))
            {
               tempMx += U[i+1][j-1][k];
               cardx++;
            }
         }
         else
         {
            tempMy += V[i][j][k];
            cardy ++;
         }
         if(flag[i][j][k-1]>= AIR_FLUID_END)
         {
            if((flag[i+1][j][k-1]>= AIR_FLUID_END)||(flag[i][j][k-1]>= AIR_FLUID_END))
            {
               tempMx += U[i+1][j][k-1];
               cardx++;
            }
            if((flag[i][j+1][k-1]>= AIR_FLUID_END)||(flag[i][j][k-1]>= AIR_FLUID_END))
            {
               tempMy += V[i][j+1][k-1];
               cardy++;
            }
         }
         tempMx /=cardx;
         tempMy /=cardy;
         deltax = tempMx - U[i][j][k];
         deltay = tempMy - V[i][j][k];
         mdeltaz = W[i][j][k] - W[i][j][k+1];
         excesszby2 = (mdeltaz - deltax - deltay)/2;
         U[i+1][j][k] = tempMx + excesszby2;
         V[i][j+1][k] = tempMy + excesszby2; 
         //cout<<"CompCell 2";
         return;
      }
   
      case 2:
      {
         if((flag[i+1][j+1][k]>= AIR_FLUID_END)|| (flag[i][j+1][k]>= AIR_FLUID_END))
         {
            tempMx += U[i+1][j+1][k];
            cardx++;
         }
         if((flag[i][j+1][k+1]>= AIR_FLUID_END)|| (flag[i][j+1][k]>= AIR_FLUID_END))
         {
            tempMz += W[i][j+1][k+1];
            cardz++;
         }
         if(flag[i-1][j][k]>= AIR_FLUID_END)
         {
            tempMx += 2* U[i][j][k] - U[i-1][j][k];
            cardx++;
            if((flag[i-1][j][k+1]>= AIR_FLUID_END)|| (flag[i-1][j][k]>= AIR_FLUID_END))
            {
               tempMz += W[i-1][j][k+1];
               cardz++;
            }
         }
         else
         {
            tempMx += U[i][j][k];
            cardx++;
         }
         if(flag[i][j-1][k]>= AIR_FLUID_END)
         {
            if((flag[i+1][j-1][k]>= AIR_FLUID_END)|| (flag[i][j-1][k] >= AIR_FLUID_END))
            {
               tempMx += U[i+1][j-1][k];
               cardx++;
            }
            if((flag[i][j-1][k+1]>= AIR_FLUID_END)|| (flag[i][j-1][k] >= AIR_FLUID_END))
            {
               tempMz += W[i][j-1][k+1];
               cardz++;
            }
         }
         if(flag[i][j][k-1]>= AIR_FLUID_END)
         {
            tempMz += 2* W[i][j][k] - W[i][j][k-1];
            cardz++;
            if((flag[i+1][j][k-1]>= AIR_FLUID_END)||(flag[i][j][k-1] >= AIR_FLUID_END))
            {
               tempMx += U[i+1][j][k-1];
               cardx++;
            }
         }
         else               
         {
            tempMz += W[i][j][k];
            cardz++;
         }            
         tempMx /=cardx;
         tempMz /=cardz;
         deltax = tempMx - U[i][j][k];
         deltaz = tempMz - W[i][j][k];
         mdeltay = V[i][j][k] - V[i][j+1][k];
         excessyby2 = (mdeltay - deltax - deltaz)/2;
         U[i+1][j][k] = tempMx + excessyby2;
         W[i][j][k+1] = tempMz + excessyby2; 
         //cout<<"CompCell 3";
         return;
      }
      case 3:
      {
         U[i+1][j][k]= -(W[i][j][k+1]+V[i][j+1][k])+(W[i][j][k]+V[i][j][k])+U[i][j][k];   
         //cout<<"CompCell 4";
         return;
      }
      case 4:
      {           
         if((flag[i+1][j+1][k]>= AIR_FLUID_END)||(flag[i+1][j][k]>= AIR_FLUID_END))
         {
            tempMy += V[i+1][j+1][k];
            cardy++;
         }
         if((flag[i+1][j][k+1]>= AIR_FLUID_END)||(flag[i+1][j][k]>= AIR_FLUID_END))
         {
            tempMz += W[i+1][j][k+1];
            cardz++;
         }
         if(flag[i][j-1][k]>= AIR_FLUID_END)
         {
            tempMy += 2* V[i][j][k] - V[i][j-1][k];
            cardy++;
            if((flag[i][j-1][k+1]>= AIR_FLUID_END)||(flag[i][j-1][k] >= AIR_FLUID_END))
            {
               tempMz += W[i][j-1][k+1];
               cardz++;
            }
         }
         else
         {
            tempMy += V[i][j][k];
            cardy++;
         }
         if(flag[i-1][j][k]>= AIR_FLUID_END)
         {
            if((flag[i-1][j+1][k]>= AIR_FLUID_END)||(flag[i-1][j][k]>= AIR_FLUID_END))
            {
               tempMy += V[i-1][j+1][k];
               cardy++;
            }
            if((flag[i-1][j][k+1]>= AIR_FLUID_END)||(flag[i-1][j][k]>= AIR_FLUID_END))
            {
               tempMz += W[i-1][j][k+1];
               cardz++;
            }
         }
         if(flag[i][j][k-1]>= AIR_FLUID_END)
         {
            tempMz += 2* W[i][j][k] - W[i][j][k-1];
            cardz++;
            if((flag[i][j+1][k-1]>= AIR_FLUID_END)||(flag[i][j][k-1]>= AIR_FLUID_END))
            {
               tempMy += V[i][j+1][k-1];
               cardy++;
            }
         }
         else
         {
            tempMz += W[i][j][k];
            cardz++;
         }
         tempMy /=cardy;
         tempMz /=cardz;
         deltay = tempMy - V[i][j][k];
         deltaz = tempMz - W[i][j][k];
         mdeltax = U[i][j][k] - U[i+1][j][k];
         excessxby2 = (mdeltax - deltay - deltaz)/2;
         V[i][j+1][k] = tempMy + excessxby2;
         W[i][j][k+1] = tempMz + excessxby2;   
         //cout<<"CompCell 5";
         return;
      }
      case 5:
      {
         V[i][j+1][k]=-(W[i][j][k+1]+U[i+1][j][k])+(W[i][j][k]+U[i][j][k])+V[i][j][k];  
         //cout<<"CompCell 6";
         return;
      }
      case 6:
     
      {
         W[i][j][k+1]=-(V[i][j+1][k]+U[i+1][j][k])+(V[i][j][k]+U[i][j][k])+W[i][j][k]; 
         //cout<<"CompCell 7";
         return;
      }
      case 7:
      {
         std::cout<<"no free faces-("<<i<<","<<j<<","<<k<<")\n";
         exit(0);
      }
   }   
}

void Navier:: NewFluidCellInitializer( int i, int j, int k)
{
   int c[6][3] ={{-1,0,0},{1,0,0},{0,-1,0},{0,1,0},{0,0,-1},{0,0,1}};
   int coord1[6]={0,0,1,1,2,2};
   int d[12][3]={{-1,-1,0},{-1,0,-1},{-1,0,1},{-1,1,0},{0,-1,-1},
                 {0,-1,1},{0,1,-1},{0,1,1},{1,-1,0},{1,0,-1},{1,0,1},{1,1,0}};
   int dc[2][12][3]={{{-1,0,0},{-1,0,0},{-1,0,0},{-1,0,0},{0,-1,0},
                 {0,-1,0},{0,1,0},{0,1,0},{1,0,0},{1,0,0},{1,0,0},{1,0,0}},                  
                 {{0,-1,0},{0,0,-1},{0,0,1},{0,1,0},{0,0,-1},
                  {0,0,1},{0,0,-1},{0,0,1},{0,-1,0},{0,0,-1},{0,0,1},{0,1,0}}};                 
   int coordzero[12]={2,1,1,2,0,0,0,0,2,1,1,2};
   int coord1dc[2][12] ={{0,0,0,0,1,1,1,1,0,0,0,0},{1,2,2,1,2,2,2,2,1,2,1,2}};
   double tempmM[3][2]={{0.0,0.0},{0.0,0.0},{0.0,0.0}};
   int cardinal = 0;
   int card[3][2];
   int coord_temp;
   bool bool_control = false;
   double const_opp[3];
   double UVW[3][3][3][3],UVWM1[3][3][3][3];
   unsigned int flag_temp = flag[i][j][k];
   double Cell_side[3]={DX[i],DY[j],DZ[k]};

   for(int l=0; l<=2;l++)
   {
      const_opp[l] = 2*(S.g[l]/S.froude)*Cell_side[l];
   }
   do
   {
      for(int l=0; l < 6;l++)
      {
         UVW[0][1+c[l][0]][1+c[l][1]][1+c[l][2]] = U[i+c[l][0]][j+c[l][1]][k+c[l][2]];
         UVW[1][1+c[l][0]][1+c[l][1]][1+c[l][2]] = V[i+c[l][0]][j+c[l][1]][k+c[l][2]];
         UVW[2][1+c[l][0]][1+c[l][1]][1+c[l][2]] = W[i+c[l][0]][j+c[l][1]][k+c[l][2]];
         UVWM1[0][1+c[l][0]][1+c[l][1]][1+c[l][2]] = U[i+c[l][0]+1][j+c[l][1]][k+c[l][2]];
         UVWM1[1][1+c[l][0]][1+c[l][1]][1+c[l][2]] = V[i+c[l][0]][j+c[l][1]+1][k+c[l][2]];
         UVWM1[2][1+c[l][0]][1+c[l][1]][1+c[l][2]] = W[i+c[l][0]][j+c[l][1]][k+c[l][2]+1];
      }
      for(int l=0; l < 6; l++)
      {
         int coord_temp = coord1[l];
         if((flag[i+c[l][0]][j+c[l][1]][k+c[l][2]] > flag_temp) &&
             (((c[l][coord_temp] < 0) && (UVWM1[coord_temp][1+c[l][0]][1+c[l][1]][1+c[l][2]]>0)) ||
             ((c[l][coord_temp] > 0) && (UVW[coord_temp][1+c[l][0]][1+c[l][1]][1+c[l][2]]<0))
               || bool_control)) 
         
         {
            cardinal ++;
            for(int m=0;m<3;m++)
            {
               if(m != coord1[l])
               {
                  tempmM[m][0] += UVW[m][1+c[l][0]][1+c[l][1]][1+c[l][2]];
                  tempmM[m][1] += UVWM1[m][1+c[l][0]][1+c[l][1]][1+c[l][2]];
               }
               else
               {
                  if(c[l][m] < 0)
                  { 
                     double UVWtemp = UVWM1[m][1+c[l][0]][1+c[l][1]][1+c[l][2]];
                     if(UVWtemp > 0)
                     {
                       double sqrtarg = UVWtemp * UVWtemp  + const_opp[m];
                       if ((sqrtarg )>0) tempmM[m][1] += sqrt(sqrtarg);
                     }
                    //tempmM[m][0] += UVWM1[m][1+c[l][0]][1+c[l][1]][1+c[l][2]];
                    //tempmM[m][1] += UVWM1[m][1+c[l][0]][1+c[l][1]][1+c[l][2]];
                  }
               else
                  { 
                     double UVWtemp =UVW[m][1+c[l][0]][1+c[l][1]][1+c[l][2]];
                     if(UVWtemp < 0)
                     {   
                        double sqrtarg = UVWtemp * UVWtemp  - const_opp[m];
                        if ((sqrtarg )>0) tempmM[m][0] -=  sqrt(sqrtarg);
                     }
                    //tempmM[m][0] += UVW[m][1+c[l][0]][1+c[l][1]][1+c[l][2]];
                    //tempmM[m][1] += UVW[m][1+c[l][0]][1+c[l][1]][1+c[l][2]];
                  }
               }
            }
         }
      }
      for(int l=0; l < 12;l++)
      {
         UVW[0][1+d[l][0]][1+d[l][1]][1+d[l][2]] = U[i+d[l][0]][j+d[l][1]][k+d[l][2]];
         UVW[1][1+d[l][0]][1+d[l][1]][1+d[l][2]] = V[i+d[l][0]][j+d[l][1]][k+d[l][2]];
         UVW[2][1+d[l][0]][1+d[l][1]][1+d[l][2]] = W[i+d[l][0]][j+d[l][1]][k+d[l][2]];
         UVWM1[0][1+d[l][0]][1+d[l][1]][1+d[l][2]] = U[i+d[l][0]+1][j+d[l][1]][k+d[l][2]];
         UVWM1[1][1+d[l][0]][1+d[l][1]][1+d[l][2]] = V[i+d[l][0]][j+d[l][1]+1][k+d[l][2]];
         UVWM1[2][1+d[l][0]][1+d[l][1]][1+d[l][2]] = W[i+d[l][0]][j+d[l][1]][k+d[l][2]+1];
      }
      if(cardinal>0)
      {
         bool_control = false;
         for(int l=0;l<=1;l++)
            for(int m=0;m<=2;m++)
               card[m][l] = cardinal;

         for(int l=0; l < 12; l++)
         {
            int coord_temp0 = coord1dc[0][l];
            int coord_temp1 = coord1dc[1][l];
            if((flag[i+d[l][0]][j+d[l][1]][k+d[l][2]]>flag_temp) &&
                ((((d[l][coord_temp0] < 0) && (UVWM1[coord_temp0][1+d[l][0]][1+d[l][1]][1+d[l][2]]>0)) ||
                ((d[l][coord_temp0] > 0) && (UVW[coord_temp0][1+d[l][0]][1+d[l][1]][1+d[l][2]]<0)))  &&
                (((d[l][coord_temp1] < 0) && (UVWM1[coord_temp1][1+d[l][0]][1+d[l][1]][1+d[l][2]]>0)) ||
                ((d[l][coord_temp1] > 0) && (UVW[coord_temp1][1+d[l][0]][1+d[l][1]][1+d[l][2]]<0)))
                || (bool_control))) 
            {
              for(int m=0;m<=1;m++)
              {
                  int coord_temp =coord1dc[m][l];
                  if ((flag[i+dc[m][l][0]][j+dc[m][l][1]][k+dc[m][l][2]]==flag_temp) &&
                     (flag[i+dc[1-m][l][0]][j+dc[1-m][l][1]][k+dc[1-m][l][2]]<=flag_temp))
                  {
                     if(d[l][coord_temp] < 0)
                     {
                        tempmM[coord_temp][0] += UVWM1[coord_temp][1+d[l][0]][1+d[l][1]][1+d[l][2]];
                        card[coord_temp][0] ++;
                     }
                     else
                     {
                        tempmM[coord_temp][1] += UVW[coord_temp][1+d[l][0]][1+d[l][1]][1+d[l][2]];
                        card[coord_temp][1] ++;
                     }
                  }
               }
            }
         }
      }
      else
      {
         std::cout<<"no fluid face neighbors-("<<i<<","<<j<<","<<k<<")\n";
         bool no_neighbor_face[3][2]={{1,1},{1,1},{1,1}};
         int  diagonal_neighbors_list[12];
         int list_count = 0;
         for(int l=0; l < 12; l++)
         {
            bool bool_temp = (flag[i+d[l][0]][j+d[l][1]][k+d[l][2]]>flag_temp);
            for (int m=coordzero[l]+1;m!=coordzero[l];m=(m+1)%3)
                 bool_temp = bool_temp && 
                (((d[l][m] < 0) && (UVWM1[m][1+d[l][0]][1+d[l][1]][1+d[l][2]]>0)) ||
                ((d[l][m] > 0) && (UVW[m][1+d[l][0]][1+d[l][1]][1+d[l][2]]<0))
                || bool_control);
            if(bool_temp)
            {
               diagonal_neighbors_list[list_count] = l;
               list_count ++;
               cardinal ++;
               coord_temp = coordzero[l];
               tempmM[coord_temp][0] += UVW[coord_temp][1+d[l][0]][1+d[l][1]][1+d[l][2]];
               no_neighbor_face[coord_temp][0] = false;
               card[coord_temp][0] ++;
               tempmM[coord_temp][1] += UVWM1[coord_temp][1+d[l][0]][1+d[l][1]][1+d[l][2]];
               no_neighbor_face[coord_temp][1] = false;
               card[coord_temp][1] ++;
               coord_temp = (coord_temp+1)%3;
               do
               {
                  if(d[l][coord_temp]== -1)
                  {
                     tempmM[coord_temp][0] += UVWM1[coord_temp][1+d[l][0]][1+d[l][1]][1+d[l][2]];
                     no_neighbor_face[coord_temp][0] = false;
                     card[coord_temp][0] ++;
                  }
                  else
                  {
                     tempmM[coord_temp][1] += UVW[coord_temp][1+d[l][0]][1+d[l][1]][1+d[l][2]];
                     no_neighbor_face[coord_temp][1] = false;
                     card[coord_temp][1] ++;
                  }              
                  coord_temp = (coord_temp + 1) %3;
               }
               while(coord_temp != coordzero[l]);
            }
         }
         for(int m=0; m < list_count; m++)
         {
            int l =  diagonal_neighbors_list[m];
            int coord_temp = (coordzero[l]+1)%3;
            do
            {
               if(d[l][coord_temp]== -1)
               {
                  if(no_neighbor_face[coord_temp][1])
                  {
                     tempmM[coord_temp][1] += UVWM1[coord_temp][1+d[l][0]][1+d[l][1]][1+d[l][2]];
                     card[coord_temp][1] ++;
                  }
               }
               else
               {  
                  if(no_neighbor_face[coord_temp][0])
                  {
                     tempmM[coord_temp][0] += UVW[coord_temp][1+d[l][0]][1+d[l][1]][1+d[l][2]];
                     card[coord_temp][0] ++;
                  }
               }              
               coord_temp = (coord_temp + 1) %3;
            }
            while(coord_temp != coordzero[l]);
         }
         if(list_count == 0)
         {
            std::cout<<"no fluid edge neighbors-("<<i<<","<<j<<","<<k<<")\n";
            int e[8][3]={{-1,-1,-1},{-1,-1,1},{-1,1,-1},{-1,1,1},{1,-1,-1},{1,-1,1},{1,1,-1},{1,1,1}};
            for(int l=0; l < 8;l++)
            {
               UVW[0][1+e[l][0]][1+e[l][1]][1+e[l][2]] = U[i+e[l][0]][j+e[l][1]][k+e[l][2]];
               UVW[1][1+e[l][0]][1+e[l][1]][1+e[l][2]] = V[i+e[l][0]][j+e[l][1]][k+e[l][2]];
               UVW[2][1+e[l][0]][1+e[l][1]][1+e[l][2]] = W[i+e[l][0]][j+e[l][1]][k+e[l][2]];
               UVWM1[0][1+e[l][0]][1+e[l][1]][1+e[l][2]] = U[i+e[l][0]+1][j+e[l][1]][k+e[l][2]];
               UVWM1[1][1+e[l][0]][1+e[l][1]][1+e[l][2]] = V[i+e[l][0]][j+e[l][1]+1][k+e[l][2]];
               UVWM1[2][1+e[l][0]][1+e[l][1]][1+e[l][2]] = W[i+e[l][0]][j+e[l][1]][k+e[l][2]+1];
            }
            for(int l=0; l < 8;l++)
            {
               bool bool_temp = (flag[i+e[l][0]][j+e[l][1]][k+e[l][2]]> flag_temp);
               for (int m=0;m<=2;m++)
                    bool_temp = bool_temp && 
                    (((e[l][m] < 0) && (UVWM1[m][1+e[l][0]][1+e[l][1]][1+e[l][2]]>0)) ||
                    ((e[l][m] > 0) && (UVW[m][1+e[l][0]][1+e[l][1]][1+e[l][2]]<0)));
               if( bool_temp || bool_control)
               {
                  diagonal_neighbors_list[list_count] = l;
                  list_count ++;
                  for(int m=0; m<=2;m++)
                  {
                     if(e[l][m]== -1)
                     {
                        tempmM[m][0] += UVWM1[m][1+e[l][0]][1+e[l][1]][1+e[l][2]];
                        no_neighbor_face[m][0] = false;
                        card[m][0] ++;
                        
                     }
                     else
                     {
                        
                        tempmM[m][1] += UVW[m][1+e[l][0]][1+e[l][1]][1+e[l][2]];
                        no_neighbor_face[m][1] = false;
                        card[m][1] ++;
                     }              
                  }           
               }
            }
            if(list_count == 0) 
            {
               std::cout<<"no fluid vertex neighbors-("<<i<<","<<j<<","<<k<<")\n"<<std::flush;
               bool_control = !(bool_control);
            }
            else
            {  bool_control = false;
               for(int n=0; n < list_count;n ++)
               {
                  int l = diagonal_neighbors_list[n];
                  for(int m=0; m<=2;m++)
                  {
                     if(e[l][m]== -1)
                     {
                        if(no_neighbor_face[coord_temp][1])
                        tempmM[m][1] += UVWM1[m][1+e[l][0]][1+e[l][1]][1+e[l][2]];
                        card[m][1] ++;
                     }
                     else
                     {
                        if(no_neighbor_face[coord_temp][0])
                           tempmM[m][0] += UVW[m][1+e[l][0]][1+e[l][1]][1+e[l][2]];
                           card[m][0] ++;
                     }              
                  }           
               }
            }
         }
         else bool_control = false;
      }
   }
   while(bool_control);
if(ISSLIP(flag[i+1][j][k])) U[i][j][k] = 0.0;
else if((flag[i-1][j][k]<= flag_temp) && (card[0][0] != 0)) U[i][j][k] = tempmM[0][0]/card[0][0];
   
if(ISSLIP(flag[i-1][j][k])) U[i+1][j][k] = 0.0;
else
{
   //if((flag[i+1][j][k]<= flag_temp) && (card[0][1]!=0)) U[i+1][j][k] = tempmM[0][1]/card[0][1];
   if((flag[i+1][j][k]< flag_temp) && (card[0][1]!=0)) U[i+1][j][k] = tempmM[0][1]/card[0][1];
   else if((flag[i+1][j][k]==flag_temp) && (card[0][1] != 0)) U[i+1][j][k] = 0.5 *((tempmM[0][1]/card[0][1])+ U[i+1][j][k]);
}

if(ISSLIP(flag[i][j+1][k])) V[i][j][k] = 0.0;
else if((flag[i][j-1][k]<=flag_temp) && (card[1][0] != 0)) V[i][j][k] = tempmM[1][0]/card[1][0];
 
if(ISSLIP(flag[i][j-1][k])) V[i][j+1][k] = 0.0;
else
{
//if((flag[i][j+1][k]<=flag_temp) && (card[1][1] != 0)) V[i][j+1][k] =tempmM[1][1]/card[1][1];
if((flag[i][j+1][k]<flag_temp) && (card[1][1] != 0)) V[i][j+1][k] =tempmM[1][1]/card[1][1];
else if((flag[i][j+1][k]==flag_temp) && (card[1][1] != 0)) V[i][j+1][k] = 0.5 *((tempmM[1][1]/card[1][1])+ V[i][j+1][k]);
}
if(ISSLIP(flag[i][j][k+1])) W[i][j][k] = 0.0;
else if((flag[i][j][k-1]<=flag_temp) && (card[2][0] != 0)) W[i][j][k] = tempmM[2][0]/card[2][0];
 
if(ISSLIP(flag[i][j][k-1])) W[i][j][k+1] = 0.0;
else
{
//if((flag[i][j][k+1]<=flag_temp) && (card[2][1] != 0)) W[i][j][k+1] = tempmM[2][1]/card[2][1];
if((flag[i][j][k+1]<flag_temp) && (card[2][1] != 0)) W[i][j][k+1] = tempmM[2][1]/card[2][1];
else if((flag[i][j][k+1]==flag_temp) && (card[2][1] != 0)) W[i][j][k+1] = 0.5 *((tempmM[2][1]/card[2][1])+ W[i][j][k+1]);
}
//cout<<"i,j,k = "<<i<<","<<j<<","<<k<<"/n";
//for(int ii=0;ii<3;ii++)
   //for(int jj=0;jj<2;jj++)
       //cout<<"Temp("<<ii<<","<<jj<<")="<<tempmM[ii][jj]<<"\n";
}




void Navier::Call_Update_Slip_Outflow_Cells()
{
   cout<< "\n Entrei na Call_Update_Slip_Outflow_Cells \n";
   int d1[6]={-1,1,0,0,0,0};
   int d2[6]={0,0,-1,1,0,0};
   int d3[6]={0,0,0,0,-1,1};
   int factor0[6]={1,1,1,1,1,1};
   int factor1[6]={0,0,0,0,0,0};
   int mn[6]={1,0,3,2,5,4};
   int i, j, k;
   cout<<"\n Slip_list \n";
   //for(int i_temp=0;i_temp<slip_list_size;i_temp++)
      //cout<<slip_list[i_temp]<<",";
   for(int i_temp=0;i_temp<slip_list_size;i_temp++)
   {
      int indice =slip_list[i_temp];
      //cout<<"\n indice"<<indice<<"\n";
      k = indice % (S.gridp[2]+2);
      j = (indice / (S.gridp[2]+2)) % (S.gridp[1]+2);
      i = indice / ((S.gridp[1]+2)*(S.gridp[2]+2)); 
      int n=(flag[i][j][k])%8;
      //cout<<"\n i,j,k,n ="<<i<<","<<j<<","<<k<<","<<n<<"\n";
      IFFLUID_INFLOW(flag[i+d1[n]][j+d2[n]][k+d3[n]]) 
      {
         int m=mn[n];
         factor0[n] = 0;
         //factor1[n] = 0;
         factor0[m] = 0;
         factor1[m] = -1;
         Update_Slip_Outflow_Cells(i+d1[n],j+d2[n],k+d3[n],i,j,k,factor0, factor1);
         /*cout<<"\n i,j,k ="<<i<<","<<j<<","<<k<<"\n";
         cout<<"\n U,V,W ="<<U[i][j][k]<<","<<V[i][j][k]<<","<<W[i][j][k]<<"\n";
         cout<<"\n U,V,WM1 ="<<U[i+1][j][k]<<","<<V[i][j+1][k]<<","<<W[i][j][k+1]<<"\n";
         cout<<"\n UM1,V,W ="<<U[1][j][k]<<","<<V[1][j][k]<<","<<W[1][j][k]<<"\n";
         cout<<"\n UM1,V,WM1 ="<<U[2][j][k]<<","<<V[1][j+1][k]<<","<<W[1][j][k+1]<<"\n";*/
         factor0[m] = 1;
         factor1[m] = 0;
         factor0[n] = 1;
         //factor1[n] = 1;
      }
   }
   for(int i_temp=0;i_temp<outflow_list_size;i_temp++)
   {
      int indice =outflow_list[i_temp];
      k = indice % S.gridp[2];
      j = (indice / S.gridp[2]) % S.gridp[1];
      i = indice / (S.gridp[1]*S.gridp[2]); 
      int n=flag[i][j][k]%8;
      IFFLUID_INFLOW(flag[i+d1[n]][j+d2[n]][k+d3[n]]) 
      {
         //int m=mn[n];
         factor0[n] = 0;
         factor1[n] = 1;
         //factor1[m] = 0;
         Update_Slip_Outflow_Cells(i+d1[n],j+d2[n],k+d3[n],i,j,k,factor0,factor1);
         factor0[n] = 1;
         factor1[n] = 0;
      }
   }
}

void Navier::Update_Slip_Outflow_Cells(int i1, int j1, int k1, int i2, int j2, int k2,int *factor0,int *factor1)
{  //(i1,j1,k1) � a posicao da celula onde est�o os dados para fazer altera��o
   //(i2,j2,k2) � a posicao da celula que ser� alterada
   U[i2][j2][k2]= factor0[0]*U[i1][j1][k1] + factor1[0]*U[i1+1][j1][k1];
   U[i2+1][j2][k2]=factor0[1]*U[i1+1][j1][k1] + factor1[1]*U[i1][j1][k1];

   V[i2][j2][k2]= factor0[2]*V[i1][j1][k1] + factor1[2]*V[i1][j1+1][k1];
   V[i2][j2+1][k2]=factor0[3]*V[i1][j1+1][k1]+ factor1[3]*V[i1][j1][k1];

   W[i2][j2][k2]= factor0[4]*W[i1][j1][k1] + factor1[4]*W[i1][j1][k1+1];
   W[i2][j2][k2+1]=factor0[5]*W[i1][j1][k1+1]+ factor1[5]*W[i1][j1][k1];
}




